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abstract 


The project, entitled "Compufcti " 
granted by NASA Langley Research Cent er ’ (N 11 .) tQ a distribute d parameter 

l^eter Nation, — suppression an, control 

££2 of large flexible aerospace structures. 

(\ ocm research was the distributed parameter modeling 
The major focus of the first-year ( ) configuration s. As an example, the Low- 

of large flexible aerospace structures wi P Satellite Model had been used as the 

power Atmospheric Compensation be summarize d as follows. A 

testbed and physical object. The main a ^P^ suppor ting frame had been built at NASA 
physical research LACE Model with y ^ ^ mo del SQ that the natural frequencies 

LaRC The classical modal test had been co analytical results obtained from 

So model were measured wWch dynamic formulation of 

the distributed parameter model of the rigid-bodv assembled structures by using transfer 

the distributed parameter model for J ^though the mathematical model must be 

matrix method had been f f " e “ffi“tures. the methodology developed m 
derived, depending on the feature of di tpther-beam-rigid-body assembled structures, 

the first-year research is definitely ” 1 “ e ° I | frequencies of the system. The analybcd 

The software PDEMOD were used to find the natu ^ report mtifled "Dtstnbuted 

values were comparable t0 1 ** ^ Model by Using Transfer Matrix Method had 
Parameter Formulation of LACE n f tbe report had been presented at the 9th 

been submitted to NASA URC,,,. A Lge SpLe Structures*,. The last two-year 

VPI&SU Symposium on Dynamics * completion of the current version 

research was conducted on two phases. The first phase ^ ^ £he dyna nuc 

of PDEMOD Code and its related documentauo . ^ P ^ ^ fj^sA Marshall Space 

properties of a two-dimensional ground-bas p a methodology for vibration 

^I^^^SparamemrmcdeUng. 

This report concentrates on die research jf j^e PDEMOD^odeTmd ^emt 

accomplishments can be summarized as fo ow . A m working on just before he died, 

completed based on several incomplete verst _ * results with those examples for 

The verification of the code had been su ® m ary of the theoretical background 

which the exact theoretical soluuons can be ob “ n “, ^reported in a technical paper 

of the package along with the « nfi “XTcSl P C»feenS ASME, Los Angeles, June 28- 
submitted to the Joint Apphed M « ha ““ ^, s MANUAL had been compiled, which mainly 
^ 0 ^ dam preparation, (2) Explanafion of the Subrouuues, (3) 

Specification of control variables. 


. wcpr two-dimensional ground-based 

Meanwhile, a theory 

manipulator f or dynamic analysis and control o axg histicate d mathematical 

S3SSS“ “““““ 

model for future modified versio 
two technical papers [4i5 ]- 


executive summary 




the finite element method ^ usually void of represent 

elements used in t e * . g that hun dreds and thousan s o _ of the computational cost 

springs. The con q acquire analytical accuracy. number of equations, there 

- 10 be ^ 

Distributed parameter mMeU^bOTg Distributed P^ter ^.taice of 

element approach for modebng ;jxg number of modal ’ { syste m 

the advantages of lhe ability to reprint the ***** “ f ^ flexib ,e 

modal order reduction, e nuations. Continuum models hav ® solar Array Flight 

dynamics in the same systcm rf ^ ^ Contro i Laboratory ( SCOLE >m- lo W . 

Xe structures, which include tte Spaecra^ station Freedom,™, and w*^ software 

- *- 

KS‘ completed for its basic functions, ‘ ^ ^ ^ NASA 

The software package ^tst 

Langley Research Center un The initial interest of the pac (he distributed 

PDEMOD package was m a f' l g enera l spacecraft configura ^ an d connected at 

model the structural dynanucs J 8 differen dal equauons is formuta wriuen in the 

prions joined at their 



t task to summarize and 

-3S5S’i 

ssi "? :**■'• •> «"» s 

tzx&rtsssstt-sss&fsi 

dements and riP d . bodl “ hi exhibit lateral bending in two ax , „„„ nM . t ed at the elements_ 

^rjss-^r % 


demS ri^ bodieS hic ^»^bending in "’^nSd - «* "' f 

A^systerrf of partial <Jiffetenti^^^^*^^ n ^™equations^f rojhton^fot^^^^^ 

aSSsSSSsS^S^^ 

mmgsss&s&z 

researchers to apply 
SmCt ”l verification of » 

(2) ExpIananon 

r«SK is^Specificationofcontioivanabies. m0( ieling technique, Dr. Taylor 

and the other researchers expec frequen cies and mode shapes, ’ d ^u CS - and (5) design 

aspects: (1) structural dam ping; (4) con ^Sand includ^i in the current 

of modal charactensUcs;^)^ ^ ^ program had been ^comp ^ is bei ng conducted, 

^rjgs^^SSSsS: 



rusticated mathematical model for 
conceived, which may provide research has been reported in two 

future modified verstons of the PDbMU 
technical papers [4>5 ]. 


enclosures 

f this report are listed as follows, which represent the research 
The enclosures of this P 

accomplishments of this project. 

, PDEM 0D code and the computed results for four verified examples, 

2 . USER'S MANUAL; demod _ A Computer Program for Distributed Parameter Esfima » 

Of nSb^er^ace Structure^, ^gi^B^y^ematics and Finable Deflecuon 

VSg£££ 'Z 1ST- «• ***"“” • 

SUGGESTION TO FUTURE RESEARCH 

As menfioned before, the PDEMOD may 

(1) — ^Eirst S 'the m'athenwfical 

^S^sjSSSS S2=*= 

SSSSSSSA. ’T— - 

important part of the whole proc^ ££ ^ by , set f PDE s. To keep 

criterion used in lhe « aTmast « 

the static stiffness were a ® J ^ distributed parameter model su , rf fic 

shucture. ^.^“‘^^rmined based on this critenon. Further, the char 

radius of gyration, etc., were tne 



♦ „ However, the stiffness of a 


n.jw - . ntaliv rnc3.i> ujLtAA 

likelihood estimator. 
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INTRODUCTION 


The software package PDEMOD was initialized by the late Lawrence W. Taylor, Jr. at 
NASA Langley Research Center during the middle of the 1980’s. The first release of his work on 
PDEMOD package was in 1987. PDEMOD was initially developed to model the structural 
dynamics of general spacecraft configurations by using the distributed parameter approach, wfoc 
insists of three-dimensional network of flexible beams and ng.d bodies. The building blocks 
from which three-dimensional configurations can be constructed consist of (1) beams, w c 
bendingin two directions, torsion, and elongation degrees of freedom, and (2) rigid bodies, wh,ch 
are connected by any network of beam elements. The full six degrees of freedom are : at towed at 
either end of the beam. Rigid bodies can be attached to the beam at any angle or body local . 
Thenmdified Bemoulli-Euler beam equation is used to represent the bending and the wave 

equatooMfor^onM^n^^^, (PDEs) . g formulated and connected at the 

elements’ boundaries based on the compatibility conditions. The equations of motion for any 
number of rigid bodies are written in the frequency domain and m terms of the coefficients oft 
sinusoidal and hyperbolic functions which comprise the mode shapes. The force and mo m 
vectors for both ends of a single beam element can be described in terms of the spatial derivatives 
of the solutions of the corresponding PDFs. Distributed parameter modds can therefore be 
generated for any three-dimensional configurations descnbable by PDEs joined at th 
boundaries. The manual labor of generating such models is therefore avoided. 

Because of Dr. Taylor’s sudden demise, it becomes an urgent task to summarize his 
research achievements and make them available to the other researchers. This work is being 
continued and this review may provide an opportunity to more researchers to .apply 
parameter modeling techniques to a variety of aerospace structures. The verification of The code 
has been conducted by comparing the results with those examples for which the exact theoretical 
solutions can be obtained, for instance, a simple beam with various boundary conditions etc. 

The theoretical derivation of the formulation and the verification of the code have been 
included in Shen, et al. m . This USER’S MANUAL concentrates on the explanation of the 
package itself. The package PDEMOD mainly includes three parts. (1) Input dat^ wluch 
specifies structural configuration, mechanical properties of the consisting beams and ngid bodies, 
and the natural frequency range one would like to search, etc., ( ) am o y o > 

which conducts the calculation specified by the formulation developed in Shen, et aJ. m . and (3) 
Subroutines necessary for completing the calculation. Since the foimulation ben t cleariy 

described in Shen, et al. m , this Manual emphasizes the first and the third parts. It should please 
the user to know that it is necessary to read Shen, et al. m , before reading this M^uak The user 
should begin by solving some of the example problems given m Shen, et al. m . The user should 

then proceed to work on their own complex problems. „ 

The USER’S MANUAL and Shen, et al. m are the basic technical specification reference 

documents for the PDEMOD package. Although these two reference documents are sufficient 
for some users, other references, e.g. [3-14], will give more details of this package^ It is o 
recommendation that users review as many of these references as possible to gam a thorough 
understanding of the PDEMOD package. 
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Neither the late Lawrence W. Taylor, Jr., nor his working colleagues the compilers of t s 
package and the authors of some of the related papers, assume any responsibility for the val y, 
accuracy, or applicability of any results. Users must verify their own results. 
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USER’S GUIDE 


A. SUBROUTINES: 

There are 17 subroutines in the package PDEMOD The subroutines . 

, • j „„ two tvnes The nine common-purpose subroutines (TYPE 

can be categorized as t yp c Tmp ttrr ady called “LODLIB.F”, which is 

qurrOTITTNES) are contained in a SUBROUTINE LioKAKx c 

the selected part of the SUBROUTINE LIBRARY “SYSPAC” (SYStems analysis progra 
PACklge) at NASA Langley Research Center (LaRC) l2| . The SYSPAC ,s a data base a LaRC 
^he Dumofe Of making experimental data and a selection of analysts algonthnts available to 

fnterestecLresearchers studying aerodynaracSj fli^t t^cha^^^iuctumMynairucs 

SYSPACand ar^prog^ed Specifically for PDEMOD. Their format is consistent with those 

of the used in the subroutines are expressed in a vector form and have a 

commo^format The first four elements of each vector are respectively: (1) the number of rows 
fit the number of columns (3) the total number of elements, and (4) the data time interval. 
( fla: Tws matrix ignition to be readily accessible in “TS ££& 

procedures, for calling numerous subroutines, in printing and in plotting. 

PDEMOD are described as follows. 


TYPE-I SUBROUTINES: 


1. SUBROUTINE ADD (P, A, Q, B, C) 

Description: Two compatible matrices (A and B) are 
respectively) and then added: C— P*A+Q B. 


multiplied by scalars (P and 


Q, 


2. SUBROUTINE MULT (A, B, C) 

Description: Multiplies two matrices. C— A B. 

3. SUBROUTINE MAKE (A, B) 

Description: Copies B in A: A=B. 

4. SUBROUTINE SET (A, n, JJ) 

Description: Creates a null matrix with II rows and JJ columns. A [OJ 


5. 

6 . 


7. 


SUBROUTINE SPIT (A, B) 

Description: Labels and lists a matrix. 

SUBROUTINE TRANS (A, B) T 

Description: Generates a matrix transpose: B=A . 

SUBROUTINE TILDA (A, B) 

Description: Forms the matrix equivalent of a cross product 


from a vector {A}3 *i, 
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0 -a 2 a 2 

B = A = a 2 0 -a l 
_-a 2 a, 0 . 

8. SUBROUTINE JUXTV (A, B, C) ....... 

Description: Combines by juxtaposition two compatible matrices m a vertical direction. 

Ta" 

c= Uj 

9. SUBROUTINE IDENT (A, II) 

Description: Forms an identity matrix A with dimension II. 


TYPE-II SUBROUTINES: 

10. SUBROUTINE AFORM (W, A, BODREF, CONFIG, L, P, Q, R, T, INRT, MASS, 
DUM, DUN, DUO, DUP) 

Description: Forms the system matrix A. (See Shen, et al.[i], Eq.34) 

11. SUBROUTINE BODFORM (CONFIG, NBEAM, BODREF, NBODY) 

Description: Determines rigid bodies’ connectivity and the reference coordinate systems. 

12. SUBROUTINE PFORM (W, L, Z, E1X, EIY, EAZ, EIP, MPL, IPL, FZ, AGKX, 
AGKY, PF, PM) 

Description: Forms the matrices Pf and Pm (see Shen, et al.[i], Eqs.16 and 17). 


13. 


SUBROUTINE QFORM (W, L, Z, EIX, EIY, EAZ, EIP, MPL, IPL, FZ, AGKX, 
AGKY, QU, QS) 

Description: Forms the matrices Q u and Q, (see Shen, et al.[i], Eqs.14 and 15). 


14. SUBROUTINE PQFORM (W, L, E, P, Q, DUM, DUN) 

Description: Combines the matrices P and Q, respectively, by juxtaposition for multi-body, 

multi-beam system. 


15. SUBROUTINE UPPER (A, DETA) 

Description: Determines the determinant of the matrix A by transforming it as a upper 
triangular matrix. 


16. SUBROUTINE DIAG (A, SHAPE) 

Description: Diagrams the mode shapes. 

17. SUBROUTINE WSEARCH (W, DW, A, DETNEW, DETOLD, WI) 

Description: Search for the values of circular natural frequencies. 
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B. VARIABLES AND ARRAYS: 

There are a number of variables and arrays are used in PDEMOD. They are defined 

below: 

VARIABLES: 

NBEAM: Number of beams. 

NBODY: Number of rigid bodies. - . 

NA=12*NBEAM, where the number 12 indicates that there are 12 unknown coefficients for eac 

beam element (see Shen, et al.[ij, Eq.13). 

Li: Length of the ith beam. . 

EEXi: Bending rigidity (EI X ) of the ith beam. 

EIYi: Bending rigidity (EI y ) of the ith beam. 

ElPi: Torsional rigidity (GI P ) of the ith beam. 

EAZi: Axial rigidity (EA Z ) of the ith beam. 

MPLi: Mass per length of the ith beam. 

IPLi: Mass moment of inertia per length of the ith beam. 

FZi: Initial tensile force for the ith tether element. 

AGKXi: Radius of gyration about x-axis of the cross-section area of the ith beam. 

AGKYi: Radius of gyration about y-axis of the cross-section area of the ith beam. 

MASSi: Mass of the ith rigid body. 

W: Circular natural frequency. 

WINC: Increment of the value of W for iteration. 

' NIW: Specified number of iteration. 


ARRAYS: 

CONFIG (NBEAM, 3): Denotes the structural configuration: Column No.l=Beam Identification 
Number (ED)’ Column No.2=Inboard Body Identification Number; Column No.3=Outboard Bo y 
Identification Number. A negative sign prior to the numbers in columns 2 and 3 indicates 
which beam is used to define body axes. For example, a two-beam, three-body system is shown 
in Figure 1. The definitions of “Inboard Body” and “Outboard Body” are depicted at the nght- 
hand sides of the body ID Numbers. The 2 by 3 matrix CONFIG for this example is numbered as 

shown in the Table 1. 


Table 1 Example of the Matrix CONFIG 


Beam’s ID Inboard Body’s ID Outboard Body’s ID 


1 

r -1 

-2 

2 

2 

-3 
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Body 3 — Outboard Body of Beam 2 



Body 1 (Earth) — Inboard Body of Beam 1 
Figure 1 Two-Beam, Three-Body System 

The negative signs prior to the numbers located at (1, 2) and (1, 3) indicate that both bodies 1 and 
2 are defined in the beam l’s coordinate system, while the negative sign prior to the number 
located at (2, 3) indicates that body 3 is defined in the beam 2’s coordinate system. 

BODREF (NBODY, 2): Defines body’s connectivity and reference coordinate. The number of 
rows of BODREF equals to the number of bodies. Each row provides the connectivity 
information of the corresponding body, consecutively. Column No.l indicates the beam’s ID, the 
connected body uses this beam’s coordinate system as the reference system. Column No.2 
indicates the mutual location between the body and the reference beam. If the body is the inboard 
body of that beam, the number equals to zero and if the body is the outboard body of that beam, 
the number equals to one. Note that BODREF is not a input array. All the numbers will be 
produced by calling SUBROUTINE BODFORM. In fact, CONFIG has provided all necessary 
information. 

Ril ( 3 , 3) : The eccentric matrix at the attachment point between the ith beam and its inboard 
body. 

RiO (3, 3): The eccentric matrix at the attachment point between the ith beam and its outboard 
body. 

R (4+18*NBEAM): The eccentric matrix containing all of matrices by juxtaposition: Ril and 
RiO, i=l, NBEAM. 

Ti (3, 3): The coordinate-transformation matrix. 

INRTi (3, 3): The mass-moment-of-inertia matrix of the ith body. 

INRT (4+9*NBODY): The mass-moment-of-inertia matrix containing all of INRTi matrices, 
i=l, NBODY, by juxtaposition. 

E (4+9*NBEAM): Beams’ mechanical-property matrix. The first four elements are used for the 
common purpose as mentioned before. Each consecutive nine numbers represent one beam’s 
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mechanical properties, respectively, that is: E(5)-EIX1, E(6)-EIY1, E(7) EAZ1, E(8) EEP1, 
E(9)=MPL1, E(10)=IPL1, E(1 1)=FZ1, E(12)=AGKX1, E(13)=AGKY1, and so on. 


MASS (4+NBODY): Body’s mass matrix. 


P [4+(3*12)*NBEAM]: Forms matrix P which consists of all the matrices P F and P M for each 
beam element consecutively by juxtaposition (see Shen, et al. m , Eqs.16 and 17). Expressed in 
matrix form, the matrix P is constructed as 


[P] = 


[Pf, Lx12 
[Pm J3X12 
[Pf,13x12 
[Pm,13x12 


For 

For 


Beaml 


Beam 2 


Q [4+(3*12)*NBEAM]: Forms matrix Q which consists of all the matrices Q„ and Q, for each 
beam element consecutively by juxtaposition (see Shen, et aLpj, Eqs.14 and 15). Expressed in 

as 


For Beaml 


For Beaml 


A [4+(12*NBEAM)*(12*NBEAM)]: Forms system matrix A which consists of all the matrices 
A f and A m for each beam element consecutively by juxtaposition, while the matrices A F and A M 
are constructed by the element matrices P F , Pm, Qu, and Q, in the way indicated in Shen, et al.jij, 
Eq.33. 

DETOLD, and DETNEW [4+(12*NBEAM)*(12*NBEAM)]: Determine the determinant of 
the system matrix A by iteration. DETOLD is the previous one, while DETNEW the up-dated 
one. When the value of the determinant is small enough to be considered as zero, the 
corresponding natural frequency is found. 


matrix form, the matrix Q is constructed 


[Q] = 


^3x12 

LxI2 

[^UjLxU 

[Q*,]3x12 


C. CONTROL VARIABLES: 

NPROB: The problem number the user chooses to solve. If NPROB=l, then PROBLEM #1 is 
solved, and so forth. The current version of PDEMOD has encluded four verification examples: 
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EXAMPLE 1 - Bending of a Cantilevered Beam; EXAMPLE 2 - Bending of a Clamped^Clamped 
Beam- EXAMPLE 3 - Bending of a Cantilevered Beam with a Tip-Mass M, EXAMPLE 4 - 
Torsion of a Cantilevered Beam. User can add his own problems into the package and assign 
corresponding problem numbers. 


IFREQ: Index for conducting natural frequency analysis. If IFREQ-1, the natural frequency of 
the system will be computed. 


ISHP: Index for conducting mode shape analysis. If ISHP=I, the mode shape functions will be 
computed. 


D. SPECIAL NUMBERS: 

The following special numbers are used in the package to define some boundary 
conditions, null mass, or infinite mass, etc., so that the tedious modal reconstruction can be 

avoided. 


1. For a null mass, or an imaginaiy mass, MASSi, at the free end of a cantilevered beam, the 
package defines MASSi=0.9*10' 8 . Correspondingly, the diagonal elements of its mass-momen - 
of-inertia matrix INRTi should be defined as the same number (0.9*10 ). 


2. For infinite mass, such as, the mass of the foundation of a cantilevered beam, MASSi, the 
nackaae defines MASSi=999999999.0 Correspondingly, the diagonal elements of its mass- 
moment-of-inertia matrix INRTi should be defined as the same number (999999999.0), 


3. For restrictions of one-direction deflection, e.g„ bending about x-axis, the package defines that 
the rigidity to resist the deflection in this direction approaches infinity, i.e., LlXi-yy^v^.u 


4. For the elements except tether element, the package defines the initial tensile force FZi 0.0 
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PDEMOD Program 


Computer-Printout Codes 


nooonno on 


123 

543 

707 

777 


PROGRAM PDEMODl 

INTEGER BODREF, CONFIG, OBODY 

REAL INRTl , INRT2 , INRT3 , INRT4 , MASS1 , MASS2 , MASS3 , MASS4 , 

1L1 , L2 , L3 , MPLl , MPL2 , MPL3 , MASS , 

9TPT.1 IPL2 I PL3 , L, INRT, INRTl, INTRO 
Dimension dum{ 13 oo), dun(1300), duo(1300) ,dup(1300) , 

2INRTl(13) , INRT2 (13) , INRT3 (13) , INRT4 (13) , 

3R1I ( 13 ) , RIO ( 13 ) , R2 1 ( 13 ) , R20 ( 13 ) , R3 1 ( 13 ) , R30 ( 13 ) , 

4E (40) , A (2 004) ,DETOLD(40) , CONFIG (5, 3) , 

5DETNEW (40), BWl ( 8 ) , BW2 ( 8 ) , BW3 ( 8 ) , 

7PF (40 ) , PM (40) , R ( 94 ) , 

9MASSU4 1 ) 'lNRT(69) BW<34) , TBODY ( 13 ) , TBEAM ( 13 ) , QO (40 ) , QS (40 ) 

raSiNSION RI(13),RO(13),INKTI(13),INRTO(13),P(436),Q(436), 

1TIT(13) ,TRT(13) , BODREF (10, 2) , Rl (13) ,CGI(7) ,CGTOT(7) , 

2TI (13) , SHAPE (2004) ,XJI (2004) ,D1 (454) , 

3DELF (40) , DELC (40) , 

4 FORMAtT2 8H 3 THE 2 SYSTEM 3 MATRIX WILL HAVE , 13 , 16HROWS AND COLUMNS) 

FORMAT (415, E15 .5) 

FORMAT (7E14 .5) 

FORMAT (7110) 

USER CAN CONTROL WHICH PROGRAM FUNCTIONS ARE ACTIVE (=1), AND NOT(=0) 

IFREQ=1 

ISHP=0 

USER MUST SELECT THE PROBLEM NUMBER DESIRED, i . e .NPROB=l , etc . 

PROBLEM #1 IS THE CANTILEVER BEAM BENDING; 

PROBLEM #2 IS THE CLAMPED-CLAMPED BEAM BENDING; 

PROBLEM #3 IS THE CANTILEVERED BEAM BENDING WITH A TIP_MASS; 

PROBLEM #4 IS THE CLAMPED-CLAMPED BEAM TORSION. 


NPROB=2 

if (nprob.eq. 1) goto 1001 
if (nprob. eq. 2 ) goto 1002 
if (nprob. eq. 3 ) goto 1003 
if (nprob. eq. 4 ) goto 1004 

c-PROBLEM 1: ********** 


c 

1001 

C 

C 

C 

C 

c 

c 

c 

C 


140 

163 


PROBLEM #1 is the CANTILEVERED BEAM BENDING (nbeam 1) . 

THE MATRIX "CONFIG" DENOTES THE STRUCTURAL CONFIGURATION 

Id! COL# 2 - INBOARD BODY I.D., COL#3=OUTBOARD BODY 
NEGATIVE SIGN DENOTES WHICH BEAM IS USED TO DEFINE BODY AXES. 

For the CANTILEVERED BEAM, BEAMl links the Bodyl (EARTH) to 
Body2 (a Null-mass) . Bodyl uses the Inboard end of Beaml to 
define its axes. Body2 uses the outboard end of Beaml to define 


its axes. 

CONFIG (1, 1) =1 
CONFIG (1, 2) =-l 
CONFIG (1, 3) =-2 
NBEAM=1 

DO 140 IBEAM = 1, NBEAM , 

PRINT 777, CONFIG (IBEAM, 1) , CONFIG (IBEAM, 2) , CONFIG (IBEAM, 3 ) 
NA=12*NBEAM 

CALL BODFORM (CONFIG, NBEAM, BODREF , NBODY) 

DO 163 I = 1, NBODY 

PRINT 777,1, BODREF (1,1), BODREF (1,2) 

CALL SET(WI, 100,1) 

WI (1) =0 . 

INPUT FOR BEAM#1, BODY#l INBOARD AND BODY#2 OUTBOARD 
Ll=130. 

E 1X1=40000000 . 

EIY1=999999999 . 0 



EIPl=999999999 . 0 
MPL1=. 09556 
IPL1= .2907 
EAZ1=999999999 . 0 
FZ1=0 . 

AGKX1=999999999 . 

AGKY1=999999999 . 

CALL SET (Rll ,3,1) 

Rll (5) =0 . 

Rll ( 6) =0 . 

Rll (7) =0 . 

CALL TILDA (Rll , DUM) 

CALL MAKE (Rll, DUM) 

CALL SET (RIO, 3,1) 

RIO (5) =0 . 0 
RIO (6) =0 . 0 
RIO ( 7 ) =0 . 0 
CALL TILDA (RIO, DUM) 

CALL MAKE (RIO, DUM) 

CALL SET (Tl ,3,3) 

T1 ( 5 ) =1 . 

Tl (9) =1 . 

Tl (13 ) =1 . 

C INPUT FOR BODY # 1 , THE FIXED END, THE EARTH. 

MASS1=999999999 . 0 
CALL SET ( INRT1 ,3,3) 

INRT1 (5) =999999999 . 0 
INRT1 (9) =999999999 . 0 
INRTl ( 13 ) =999999999.0 
PRINT 707 , MASS1 

C INPUT FOR BODY# 2, THE NULL -MASS. 

MASS2=0. 00000009 
CALL SET ( INRT2 ,3,3) 

INRT2 (5) =0.00000009 
INRT2 (9) =0. 00000009 
INRT2 (13) =0.00000009 
PRINT 7 07, MAS S2 
CALL SET ( E , NBEAM , 9 ) 

E (5) =EIXl 

E ( 6) =EIYl 

E ( 7 ) =EAZl 

E ( 8 ) =EIP1 

E (9) =MPL1 

E ( 10 ) =IPLl 

E ( 11) =FZ1 

E ( 12 ) =AGKX1 

E ( 13 ) =AGKY1 

CALL SET (L, NBEAM, 1) 

L (5) =L1 

CALL SET (MASS, NBODY, 1) 

MASS ( 5 ) =MASS1 
MASS ( 6 ) =MASS2 
CALL MAKE (INRT, INRTl) 

CALL JUXTV ( INRT, INRT2, INRT) 

CALL MAKE (T, Tl) 

CALL MAKE (R, Rll) 

CALL JUXTV (R, RIO, R) 

W=4 . 0 
WINC=1 . 01 
NIW=600 
GO TO 1000 

•PROBLEM 2: ********** 
c 

1002 CONTINUE 

C PROBLEM #2 is the CLAMPED_CLAMPED BEAM BENDING (nbeam=l) . 



C THE MATRIX "CONFIG" DENOTES THE STRUCTURAL CONFIGURATION 

C COL# 1 =BEAM I.D., COL#2 -INBOARD BODY I.D., COL#3=OUTBOARD BODY 

C NEGATIVE SIGN DENOTES WHICH BEAM IS USED TO DEFINE BODY AXES, 

c For the CANTILEVERED BEAM, BEAMl links the Bodyl (EARTH) to 

c Body2 (EARTH, EITHER) . Body! uses the Inboard end of Beaml to 

- define its axes. Body2 uses the outboard end of Beaml to define 

its axes. 

CONFIG (1, 1 ) =1 
CONFIG (1, 2) =-l 
CONFIG (1, 3 )=-2 
NBEAM=1 

DO 240 IBEAM = 1 , NBEAM 

240 PRINT 777, CONFIG ( IBEAM, 1 ), CONFIG ( IBEAM, 2 ), CONFIG ( IBEAM, 3 ) 

NA=12 * NBEAM 

CALL BODFORM (CONFIG, NBEAM, BODREF , NBODY) 

DO 263 I = 1, NBODY 

263 PRINT 777,1, BODREF (1,1) , BODREF (1,2) 

CALL SET(WI, 100, 1) 

WI(1)=0. 

C INPUT FOR BEAM# 1 , BODY#l INBOARD AND BODY#2 OUTBOARD 

Ll=130 . 

EIX1=40000000 . 

EIY1=99 9999999 . 0 
EIP1=999999999 . 0 
MPL1=. 09556 
IPL1= . 2907 
EAZ1=999999999 . 0 
FZ1=0 . 

AGKXl=9 999 . 

AGKY1=9 99999999 . 

CALL SET (Rll , 3,1) 

Rll (5 ) =0 . 

Rll (6) =0 . 

Rll (7 ) =0 . 

CALL TILDA (Rll, DUM) 

CALL MAKE (Rll, DUM) 

CALL SET (RIO, 3,1) 

RIO (5) =0 . 0 
RIO ( 6 ) =0 . 0 
RIO (7 ) =0 . 0 
CALL TILDA (RIO, DUM) 

CALL MAKE (RIO, DUM) 

CALL SET (T1 , 3 , 3 ) 

T1 (5) =1 . 

T1 (9 ) =1 . 

T1 (13 ) =1 . 

C INPUT FOR BODY # 1 , THE FIXED END, THE EARTH. 

MASS1=999999999 . 0 
CALL SET ( INRTl ,3,3) 

INRT1 (5) =9999999.0 
INRTl (9) =9999999.0 
INRTl (13) =9999999.0 
PRINT 707 , MASS1 

C INPUT FOR BODY# 2 , THE NULL -MAS S . 

MASS2=999999999 . 0 
CALL SET ( INRT2 ,3,3) 

INRT2 (5) =9999999.0 
INRT2 (9) =9999999.0 
INRT2 (13) =9999999.0 
PRINT 707 , MASS2 
CALL SET (E, NBEAM, 9) 

E ( 5 ) =EIXl 
E (6) =EIY1 
E (7 ) =EAZl 
E ( 8 ) =EIPl 
E ( 9 ) =MPL1 



E ( 10 ) =IPLl 

E(11)=FZ1 

E ( 12 ) =AGKXl 

E ( 13 ) =AGKYl 

CALL SET ( L , NBEAM , 1 ) 

L (5) =Ll 

CALL SET (MASS , NBODY , 1 ) 

MASS (5) =MASS1 
MASS ( 6 ) =MASS2 
CALL MAKE ( INRT , INRTl ) 

CALL JUXTV ( INRT , INRT2 , INRT ) 
CALL MAKE (T, Tl) 

CALL MAKE (R, R1I) 

CALL JUXTV (R, RIO, R) 

W=25 . 0 
WINC=1 . 01 
NIW=300 
GO TO 1000 

c-PROBLEM 3 ********** 


c 

1003 

C 

c 

C 

C 

C 

c 

c 

c 

C 


340 

363 

C 


PROBLEM E #3 is the CANTILEVERED BEAM BENDING (nbeam=l) with a 

DENOTES THE STRUCTURAL CONFIGURATION 
COLftl^BEAM I.D , COL# 2 -INBOARD BODY I.D., COL#3=OUTBOARD BODY 
NEGATIVE SIGN DENOTES WHICH BEAM IS USED TO DEFINE BODY AXES. 

For the CANTILEVERED BEAM, BEAMl links the Bodyl (EARTH) to 
Bodv2 Bodvl uses the Inboard end of Beaml to . 

define i?s axes. Body2 uses the outboard end of Beaml to define 


its axes . 
CONFIG (1, 1) =1 
CONFIG (1, 2) =-l 
CONFIG (1, 3) =-2 


NBEAM=1 

PRINT°777, CONFIG (IBEAM, 1) , CONFIG ( IBEAM, 2 ) , CONFIG ( IBEAM, 3 ) 


NA=12*NBEAM mrtTW . 

CALL BODFORM ( CONFIG , NBEAM , BODREF , NBODY ) 

DO 363 I = 1, NBODY 

PRINT 777,1, BODREF (1,1), BODREF (1,2) 

CALL SET (WI, 100,1) 

INPUT FOR BEAM#1 , BODY#l INBOARD AND BODY#2 OUTBOARD 


Ll=3 . 077 
EIX1=175 .9644 
EIYl=999999999 . 0 
EIP1=999999999 . 0 
MPL1=0 . 012 
IPL1=1 . OelO 
EAZl=999999999 . 0 


FZ1=0 . 

AGKXl=999999999 . 
AGKY1=999999999 . 
CALL SET(R1I, 3 , 1) 
Rll (5) =0 . 

Rll (6) =0 . 

Rll (7 ) =0 . 

CALL TILDA (Rll, DUM) 
CALL MAKE (Rll, DUM) 
CALL SET (RIO, 3,1) 
RIO (5) =0 . 0 
RIO (6) =0 . 0 
RlO (7 ) =-0.07 
CALL TILDA (RIO, DUM) 
CALL MAKE (RIO, DUM) 



c 


CALL SET (Tl ,3,3) 

T1 (5) =1 . 

Tl (9) =1. 

Tl (13) =1 . 

INPUT FOR BODY#l , THE CLAMPED_END , THE EARTH. 

MASS1=99 9999999 . 0 
CALL SET ( INRT1 ,3,3) 

INRT1 (5) =999999999 . 0 
INRTl (9) =999999999 . 0 
INRT1 (13 ) =999999999.0 
PRINT 707 ,MASS1 

INPUT FOR BODY # 2 , THE REFLECTOR. 

MASS2=4. 952/32. 2 
CALL SET ( INRT2 ,3,3) 

INRT2 (5)=2.341e-4 
INRT2 (9 ) =2 . 341e-4 
INRT2 (13 ) =1 . 524e-3 
PRINT 707 ,MASS2 
CALL SET ( E , NBEAM , 9 ) 

E (5 ) =EIXl 
E (6) =EIYl 
E (7 ) =EAZl 
E ( 8 ) =EIPl 
E (9 ) =MPL1 
E ( 10 ) =IPLl 
E (11) =FZl 
E (12) =AGKX1 
E ( 13 ) =AGKY1 
CALL SET (L, NBEAM, 1) 

L (5 ) =L1 

CALL SET (MASS, NBODY, 1) 

MASS (5) =MASS1 
MASS ( 6 ) =MASS2 
CALL MAKE ( INRT , INRTl ) 

CALL JUXTV ( INRT , INRT2 , INRT ) 

CALL MAKE (T, Tl ) 

CALL MAKE (R, Rll ) 

CALL JUXTV (R, RIO, R) 

W=10.0 
winc=l . 01 
niw=150 
GO TO 1000 

c-PROBLEM 4 ********** 
c 

1004 CONTINUE , , 

C PROBLEM #4 is the TORSION Of A CLAMOED- CLAMPED BEAM (nbeam=l) 

c WITHOUT BODY CONNECTED. 

C THE MATRIX "CONFIG" DENOTES THE STRUCTURAL CONFIGURATION 

C COL#l=BEAM I.D., COL#2 -INBOARD BODY I.D., COL#3=OUTBOARD BODY 

C NEGATIVE SIGN DENOTES WHICH BEAM IS USED TO DEFINE BODY AXES, 

c For the CANTILEVERED BEAM, BEAMl links the Bodyl (EARTH) to 

c Body2 (another clamped end) . Bodyl uses the Inboard end of Beaml to 

c define its axes. Body2 uses the outboard end of Beaml to define 

C its axes. 

CONFIG (1, 1) =1 
CONFIG (1, 2) =-l 
CONFIG (1, 3 )=-2 
NBEAM=1 

DO 440 IBEAM = 1 , NBEAM 

440 PRINT 777, CONFIG (IBEAM, 1) , CONFIG (IBEAM, 2) , CONFIG (IBEAM, 3) 

NA=12 *NBEAM 

CALL BODFORM ( CONFIG , NBEAM , BODREF , NBODY ) 
do 463 i=l,nbody 

print 777, i , bodref ( i , 1 ) , bodref ( i , 2 ) 

CALL SET (WI, 100,1) 


463 



WI (1) =0 . 

c INPUT FOR BEAM 1, BODY#l INBOARD, BODY#2 OUTBOARD 

Ll=130. 

EIX1=999999999 . 0 
EIyl=999999999 . 0 
EIpl=400000000 . 0 
MPL1=0 .09556 
IPL1= .2907 
EAZ1=999999999 . 0 
FZ1=0. 

AGKX1=9 9999999 9 . 0 
AGKyl=999999999 . 0 
CALL SET (Rll ,3,1) 

Rll (5) =0 . 

Rll (6) =0 . 

Rll (7 ) =0 . 

CALL TILDA (Rll, DUM) 

CALL MAKE (Rll, DUM) 

CALL SET (RIO, 3,1) 

RIO (5) =0 . 

RlO(6) =0 . 

RlO (7 ) =0 . 

CALL TILDA (RIO, DUM) 

CALL MAKE (RIO, DUM) 

CALL SET (Tl ,3,3) 

Tl (5) =1 . 

Tl (9) =1 . 

Tl (13 ) =1 . 

C INPUT FOR BODY l(one clamped end) 

MASS1=999999999 . 

CALL SET ( INRTl ,3,3) 

INRT1 (5) =999999999 . 0 
INRTl (9) =999999999 . 0 
INRTl (13) =999999999 . 0 
PRINT 707, MASS1 

C INPUT FOR BODY 2 (another clamped end) 

MASS2=999999999 . 0 
CALL SET ( INRT2 ,3,3) 

INRT2 (5) =999999999 . 0 
INRT2 (9) =999999999 . 0 
INRT2 (13 ) =999999999 . 0 
PRINT 707, MASS2 
CALL SET ( E , NBEAM , 9 ) 

E ( 5 ) =EIX1 

E ( 6 ) =EIY1 

E (7 ) =EAZ1 

E (8) =EIP1 

E (9) =MPL1 

E ( 10 ) =IPL1 

E (11) =FZl 

E (12) =AGKX1 

E ( 13 ) =AGKY1 

CALL SET (L, NBEAM, 1) 

L ( 5 ) =Ll 

CALL SET (MASS , NBODY , 1 ) 

MASS (5) =MASS1 
MASS ( 6 ) =MASS2 
CALL MAKE (INRT, INRTl) 

CALL JUXTV ( INRT , INRT2 , INRT ) 

CALL MAKE (T, Tl ) 

CALL MAKE (R, Rll) 

CALL JUXTV (R, RIO, R) 

W=280 . 0 
WINC=1 . 01 
NIW=300 
GO TO 1000 



c 

c ***** the common part of the main program ***** 

c 

1000 CALL SPIT (L, 2H L) 

CALL SPIT (E, 2H E) 

CALL SPIT (R, 2H R) 

CALL SPIT (T, 2H T) 

CALL SPIT (MASS, 5H MASS) 

CALL SPIT ( INRT, 5H INRT) 

CALL SET { DETNEW , NA , 1 ) 

DO 1 IW=1 , NIW 
DW=W* (WINC-1 . ) 

W=W+DW 

C FORM FORCE AND MOMENT MATRICES "PF" AND "PM" 

C FORM LINEAR AND ANGULAR DEFLECTION MATRICES "QU" AND "QS" 

CALL PQFORM (W, L, E, P, Q, DUM, DUN) 

CALL SET (A, NA, NA) 

CALL AFORM(W, A, BODREF, CONFIG, L, P, Q, R, T, INRT, MASS, 

1DUM, DUN, DUO , DUP ) 

CALL ADD ( 1 . , A, 0 . , A, A) 

CALL WSEARCH ( W , DW , A , DETNEW , DETOLD , WI ) 

1 CONTINUE 

AD=l./6. 283185 

CALL ADD (AD, WI , 0 . , WI , FI ) 

CALL SPIT (FI, 3H FI) 

IF(ISHP) 998,998,31 
3 1 CONTINUE 

NW=WI ( 1 ) + . 0 1 
DO 3 IW=1 , NW 
W=WI (IW+4) 

CALL PQFORM(W, L, E, P, Q, DUM, DUN) 

CALL AFORM(W,A, BODREF, CONFIG, L, P,Q,R,T, INRT, MASS, 

1DUM, DUN, DUO, DUP) 


92 

CALL UPPER (A, DETNEW) 
DETN=1 . 

DO 92 I = 1,NA 
DETN=DETN* DETNEW ( 1+4 ) 



CALL SET (DUM, 1 , NA) 
DUM (NA+4 ) =1 . 

CALL DIAG (A, DUM) 


30 

IF ( IW- 1)29,30,29 
CALL MAKE (SHAPE, DUM) 


29 

GO TO 28 

CALL JUXTV( SHAPE, DUM, 

SHAPE) 

28 

CONTINUE 


3 

CONTINUE 


999 

IF ( ISHP )998,998,999 
CALL SPIT ( SHAPE , 5HSHAPE ) 

998 

STOP 


c 

END 

★★★★★ 

* * * * * 

c 

SUBROUTINES 


c 

★ * ★ * ★ ★ 

★ * * * * 


SUBROUTINE AFORM (W, A, BODREF , CONFIG , L, P, Q, R, T, INRT , MASS , 
1DUM, DUN, DUO, DUP) 

INTEGER APPEND, BEAMl , BODREF, CONFIG, OBODY 
REAL INRT, INRTI, INRTO, MASS, L 

DIMENSION DUM (13 00) , DUN (13 00) , DUO (1300) , DUP (1300) , 

2 INRT (69) , INRTI (13) , MASS (14) , L (9) , R(94) ,RI (13) ,R1 (13) , 
3A(2004) , P (436) , PF (40 ) , PM (40 ) , Q (43 6 ) , QU (40) , QS (40) , 
4BODREF (10,2) , CONFIG (5, 3) ,T(49) , TBEAM ( 13 ) ,T1(13) , 

5TIT ( 13 ) , TRT ( 13 ) ,QUl(40) ,QSl(40) , TI (13 ) , RO (13 ) , 

6 INRTO (13 ) 

707 FORMAT (7E15. 4) 

NBODY=MASS ( 1 ) + . 01 



NBEAM=L ( 1 ) + . 01 

NA=NBEAM*12 

APPEND=NBODY 

DO 2 IBEAM = 1 , NBEAM 

CALL SET (TBEAM, 3,3) 

CALL SET (RI ,3,3) 

CALL SET (RO, 3,3) 

DO 13 IBLOCK = 1,9 

TBEAM ( IBLOCK+4 ) =T ( (IBEAM*3-3 ) *3+IBLOCK+4) 

RI (IBLOCK+4) =R( (IBEAM*6-6) *3+IBLOCK+4) 

13 RO (IBLOCK+4) =R( ( IBEAM*6-3 ) *3+IBLOCK+4 ) 

C FORM MATRIX "A" FOR INBOARD BODY 

CALL SET (QU, 3, 12) 

CALL SET(QS,3, 12) 

DO 12 IBLOCK = 1,36 

QU ( IBLOCK+4 ) =Q( ( IBEAM*12-12 ) *12+IBLOCK+4 ) 

12 QS( IBLOCK+4 ) =Q ( ( IBEAM*12 -9 ) *12+IBLOCK+4 ) 

CALL TRANS (TBEAM, TIT) 

IBODY=CONFIG ( IBEAM, 2 ) 

IF (IBODY) 6, 4, 4 

6 IBODY=- IBODY 
ENREF=1 . 

GO TO 18 
4 ENREF=0 . 

BEAMl =BODREF ( IBODY, 1 ) 

IO=BODREF ( IBODY, 2 ) 

IF (BEAMl -IBEAM) 70,71,70 

70 CALL SET (Tl ,3,3) 

CALL SET (QUl ,3,12) 

CALL SET (QS1 , 3,12) 

CALL SET (Rl ,3,3) 

DO 72 IBLOCK = 1,9 

Tl (IBLOCK+4) =T( (BEAMl *3 -3) *3+IBLOCK+4) 

-'■> Rl (IBLOCK+4) =R( (BEAMl*6 + 3 *10-6 ) *3 + IBLOCK+4) 

DO 73 IBLOCK = 1,36 

QUl (IBLOCK+4) =Q( (BEAM1*12+6*I0-12 ) *12+IBLOCK+4) 

73 QS1 (IBLOCK+4) =Q( (BEAMl*12+6*IO-9 ) *12+IBLOCK+4 ) 
CALL MULT (Tl, QUl, DUN) 

CALL ADD ( 1 . , RI , -1 . , Rl , DUM) 

CALL MULT ( DUM, Tl, DUO) 

CALL MULT (DUO, QS1, DUM) 

CALL ADD ( -1 . , DUN, -1 . , DUM, DUO) 

CALL MULT (TIT, DUO, DUP) 

CALL MULT (T1,QS1, DUM) 

CALL MULT (TIT, DUM, DUN) 

CALL ADD ( -1 . , DUN, 0 . , DUN, DUN) 

CALL JUXTV ( DUP , DUN , DUO ) 

CALL JUXTV ( QU , QS , DUP ) 

DO 75 IBLOCK =1,6 
DO 74 JBLOCK = 1,12 

ID= (APPEND* 6+ IBLOCK- 1) *NA+JBLOCK+BEAMl*12-8 

74 A ( ID) =DUO { (IBLOCK-1) *12+ JBLOCK+4 ) 

DO 75 JBLOCK =1,12 

ID= ( APPEND* 6 + IBLOCK- 1 ) *NA+ JBLOCK+ IBEAM* 12 - 8 

7 5 A ( ID) =DUP ( (IBLOCK-1) *12+ JBLOCK+4 ) 

APPEND=APPEND+1 

7 1 CONTINUE 

18 CONTINUE 

EM=MASS ( IBODY+4 ) 

CALL SET(INRTI, 3,3) 

DO 14 IBLOCK = 1,9 

■■ 4 INRTI (IBLOCK+4) =INRT( (IBODY*3-3 ) *3 + IBLOCK+4) 

CALL SET (PF, 3 , 12 ) 

CALL SET ( PM, 3,12) 

DO 31 IBLOCK = 1,36 

PF ( IBLOCK+4 ) =P( (IBEAM*12-12) *12+ IBLOCK+4) 



31 PM ( IBLOCK+4 ) =P ( (IBEAM* 12 -9) *12+IBLOCK+4) 

W2=W*W 
W2IN=1. /W2 

CALL MULT (RI , TBEAM, DUM) 

CALL MULT (TIT, DUM, TRT) 

FIRST BLOCK, FORCE EQUATIONS, INBOARD 
IF (EM-999999999 . ) 86, 87, 86 

87 AD=ENREF 
AE=0. 

GO TO 88 

86 AD=ENREF * EM 

AE=W2IN 

88 CONTINUE 

CALL ADD (AD, QU , AE, PF , DUO) 

CALL MULT (TRT, QS, DUN) 

CALL ADD ( 1 . , DUO , + AD , DUN , DUO ) 

C SECOND BLOCK, MOMENT EQUATIONS, INBOARD 

CALL MULT (RI, TBEAM, DUM) 

CALL MULT (DUM, PF, DUN) 

CALL MULT (TBEAM, PM, DUM) 

CALL ADD ( +W2 IN , DUM , +W2 IN , DUN , DUP ) 

CALL MULT (TBEAM, QS, DUM) 

CALL MULT (INRTI, DUM, DUN) 

CALL ADD ( 1 . , DUP , ENREF , DUN , DUP ) 

CALL JUXTV ( DUO , DUP , DUO ) 

DO 5 IBLOCK =1,6 
DO 5 JBLOCK =1,12 

5 A ( ( IBODY* 6 + IBLOCK-7 ) *NA+IBEAM*12+ JBLOCK-8 ) =DUO ( (IBLOCK-1) 

1*12 + JBLOCK+ 4 ) 

C FORM MATRIX "A" FOR OUTBOARD BODY 

CALL SET ( QU , 3,12) 

CALL SET (QS, 3 , 12 ) 

DO 19 IBLOCK = 1,36 

QU ( IBLOCK+4 ) =Q( (IBEAM*12-6) *12+ IBLOCK+4) 

, QS( IBLOCK+4 )=Q( (IBEAM*12-3 ) *12+IBLOCK+4) 

OBODY=CONFIG ( IBEAM, 3 ) 

IF (OBODY) 15,15, 16 

15 ENREF =1 . 

OBODY=- OBODY 
GO TO 17 

16 ENREF =0 . 

BEAMl=BODREF (OBODY, 1 ) 

IO=BODREF (OBODY, 2 ) 

IF (BEAM1- IBEAM) 80, 81, 80 
80 CALL SET (Tl ,3,3) 

CALL SET(QU1, 3, 12) 

CALL SET (QS1 , 3,12) 

CALL SET (Rl ,3,3) 

DO 82 IBLOCK =1,9 

Tl (IBLOCK+4) =T( ( BEAMl * 3 - 3 ) * 3 + IBLOCK+ 4 ) 

82 Rl (IBLOCK+4) =R( (BEAMl*3-3 ) *3+IBLOCK+4) 

DO 83 IBLOCK = 1,36 

QU1 (IBLOCK+4) =Q( (BEAM1*12 + 6*I0-12 ) *12+IBLOCK+4) 

83 QS1( IBLOCK+4 ) =Q ( (BEAMl* 12+6* IO-9 ) *12+IBLOCK+4) 

CALL MULT ( Tl , QUl , DUN ) 

CALL ADD ( 1 . , RI , -1 . , Rl , DUM) 

CALL MULT (DUM, Tl, DUO) 

CALL MULT (DUO, QS1, DUM) 

CALL ADD ( - 1 . , DUN , - 1 . , DUM , DUO ) 

CALL MULT (TIT, DUO, DUP) 

CALL MULT (T1,QS1, DUM) 

CALL MULT (TIT, DUM, DUN) 

CALL ADD ( - 1 . , DUN , 0 . , DUN , DUN ) 

CALL JUXTV (DUP, DUN, DUO) 

CALL JUXTV (QU,QS, DUP) 

DO 85 IBLOCK = 1,6 
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DO 84 JBLOCK =1,12 

ID= (APPEND* 6 +IBLOCK-1) *NA+ JBLOCK+BEAMl*12-8 
84 A ( ID) =DUO ( (IBLOCK-1) *12+ JBLOCK+4 ) 

DO 85 JBLOCK =1,12 

ID= ( APPEND* 6 + IBLOCK- 1 ) *NA+ JBLOCK+ IBEAM* 12 - 8 
~5 A (ID) =DUP ( (IBLOCK-1) * 12+ JBLOCK+4 ) 

APPEND=APPEND+ 1 
81 CONTINUE 

17 CONTINUE 

EM=MAS S ( OBOD Y + 4 ) 

CALL SET ( INRTO ,3,3) 

DO 24 IBLOCK =1,9 

24 INRTO ( IBLOCK+4 ) =INRT( (OBODY*3-3 ) *3+IBLOCK+4) 

CALL SET (PF, 3, 12) 

CALL SET ( PM, 3,12) 

DO 33 IBLOCK = 1,36 

PF( IBLOCK+4 )=P( (IBEAM*12-6) *12+IBLOCK+4) 

33 PM ( IBLOCK+4 ) =P ( (IBEAM*12-3) *12+IBLOCK+4 ) 

W2=W*W 

W2IN=1./W2 

CALL MULT ( RO , TBEAM , DUM ) 

CALL MULT (TIT, DUM, TRT) 

C FIRST BLOCK, FORCE EQUATIONS, OUTBOARD 

IF (EM-999999999 .) 96, 97, 96 

97 AD=ENREF 
AE=0 . 

GO TO 98 

96 AD= ENREF * EM 

AE=W2 IN 

98 CONTINUE 

CALL ADD (AD, QU, AE, PF, DUO) 

CALL MULT (TRT, QS, DUN) 

CALL ADD ( 1 . , DUO , +AD , DUN , DUO ) 

SECOND BLOCK, MOMENT EQUATIONS, OUTBOARD 
CALL MULT (RO, TBEAM, DUM) 

CALL MULT (DUM, PF, DUN) 

CALL MULT (TBEAM, PM, DUM) 

CALL ADD ( +W2 IN , DUM , +W2 IN , DUN , DUP ) 

CALL MULT (TBEAM, QS, DUM) 

CALL MULT (INRTO, DUM, DUN) 

CALL ADD ( 1. , DUP, ENREF, DUN, DUP) 

CALL JUXTV ( DUO , DUP , DUO ) 

DO 7 IBLOCK =1,6 
DO 7 JBLOCK =1,12 

7 A ( (OBODY*6+IBLOCK-7 ) *NA+ IBEAM* 12 + JBLOCK- 8 ) =DUO ( (IBLOCK-1) 

1*12+ JBLOCK+4) 

2 CONTINUE 

RETURN 
END 
c 

SUBROUTINE BODFORM (CONFIG, NBEAM, BODREF , NBODY) 

INTEGER CONFIG, BODREF 

DIMENSION CONFIG (5,3), BODREF (10,2) 

777 FORMAT (7110) 

JMAX=1 

DO 50 IBEAM = 1, NBEAM 
J=CONFIG ( IBEAM, 2 ) 

J2=J*J 

IF ( JMAX* JMAX- J2 ) 1 , 2 , 2 

1 JMAX2=J2 

2 CONTINUE 

IF ( J) 51,51, 52 
J=-J 

BODREF ( J, 1) = IBEAM 
BODREF (J, 2) =0 
52 J=CONFIG ( IBEAM , 3 ) 



J2=J*J 

IF ( JMAX*JMAX-J2 ) 3 , 4 , 4 

3 JMAX2=J2 

4 CONTINUE 

IF ( J) 53 , 53 , 54 
^3 J=-J 

BODREF ( J, 1) =IBEAM 
BODREF ( J , 2 ) = 1 
54 CONTINUE 

50 CONTINUE 

AD=JMAX2 

NBODY=SQRT (AD) + . 01 
DO 63 I = 1 , NBODY 

63 PRINT 777,1, BODREF (1,1), BODREF (1,2) 

RETURN 

END 

SUBROUTINE PFORM(W, L, Z, EIX, EIY, EAZ , EIP, MPL, IPL, FZ , AGKX, AGKY, 
1PF, PM) 

REAL L, MPL, IPL 
DIMENSION PF (40), PM(40) 

707 FORMAT (7E14. 5) 

W2=W*W 
W2IN=1. /W2 

IF (EIX-FZ*FZ* .071)8,8,9 

8 ARG=MPL/ (FZ-EIX*MPL*W2 /AGKX) 

BXAB=W* SQRT (ARG) 

BXCD=18 . 42 /L 
GO TO 10 

9 CONTINUE 

AD= . 5* (MPL*W2 /AGKX-FZ/EIX) 

AE=MPL*W2/EIX 
ARG=AD+SQRT (AD*AD+AE ) 

BXAB= SQRT (ARG) 

ARG =-AD+ SQRT (AD*AD+AE) 

BXCD= SQRT (ARG) 

10 CONTINUE 

IF (EIY-FZ*FZ* . 071)11, 11, 12 

11 ARG=MPL/ (FZ-EIY*MPL*W2 /AGKY) 

BYAB=W* SQRT ( ARG ) 

BYCD=18 . 42/L 
GO TO 13 

12 CONTINUE 

AD= . 5* (MPL*W2 /AGKY-FZ/EIY) 

AE=MPL*W2 /EIY 
ARG=AD+SQRT (AD*AD+AE ) 

BYAB= SQRT (ARG) 

ARG= - AD+ SQRT (AD*AD+AE) 

BYCD= SQRT (ARG) 

13 CONTINUE 
BZ=SQRT (MPL /EAZ ) *W 
BP=SQRT ( IPL/EIP) *W 

C FORCE MATRIX FOR BEAM 

CALL SET (PF, 3,12) 

BXAB2 =BXAB*BXAB 
BXCD2=BXCD*BXCD 
BXAB3 =BXAB2 * BXAB 
BXCD3=BXCD2*BXCD 
BYAB2 =BYAB*BYAB 
BYCD2 =B YCD *B YCD 
BYAB3 =BYAB2 *BYAB 
BYCD3 =BYCD2 *BYCD 
IF (Z) 2 , 3 , 2 

J PF ( 5 ) =+BXAB3 *EIX+BXAB*FZ 

PF (6) =0 . 

PF (7 ) =-BXCD3 *EIX-BXCD*FZ 
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PF (8) =0. 

PF ( 2 1 ) =+BYAB3*EIY+BYAB*FZ 
PF ( 2 2 ) = 0 . 

PF (23) =-BYCD3*EIY-BYCD*FZ 
PF ( 2 4 ) = 0 . 

PF (37 ) =BZ*EAZ 
PF ( 3 8 ) = 0 . 

GO TO 4 

PF ( 5 ) =-BXAB3 *EIX*COS (BXAB*L) -BXAB*FZ*COS (BXAB*L) 

PF ( 6 ) = BXAB3 *EIX*SIN (BXAB*L) +BXAB*FZ*SIN (BXAB*L) 

PF (7 ) = BXCD3*EIX*COSH(BXCD*L) -BXCD*FZ*COSH (BXCD*L) 

PF ( 8 ) = BXCD3 *EIX*SINH (BXCD*L) -BXCD*FZ*SINH (BXCD*L) 

PF (21) =-BYAB3*EIY*COS (BYAB*L) -BYAB*FZ*COS (BYAB*L) 

PF (22 ) = BYAB3 *EIY*SIN (BYAB*L) +BYAB*FZ*SIN (BYAB*L) 

PF ( 2 3 ) = BYCD3 *EI Y*COSH ( BYCD* L ) -BYCD*FZ *COSH ( BYCD*L) 

PF ( 2 4 ) = BYCD3 *EIY*SINH (BYCD*L) -BYCD*FZ*SINH (BYCD*L) 

PF (37) =-BZ*EAZ*COS (BZ*L) 

PF (38) =BZ*EAZ*SIN (BZ*L) 

4 CONTINUE 

C MOMENT MATRIX FOR BEAM 

CALL SET (PM, 3, 12) 

IF (Z) 5, 6, 5 

6 PM (17) =0 . 

PM ( 18 ) =+BXAB2 *EIX+FZ 
PM (19) =0 . 

PM (20) =-BXCD2 *EIX+FZ 
PM ( 9 ) = 0 . 

PM (10) =-BYAB2*EIY-FZ 
PM (11) =0 . 

PM ( 12 ) =+BYCD2 *EIY-FZ 
PM(39) =BP*EIP 
PM ( 4 0 ) = 0 . 

GO TO 7 

PM(17) =-BXAB2*EIX*SIN(BXAB*L) -FZ*SIN (BXAB*L) 

PM (18) =-BXAB2*EIX*COS (BXAB*L) -FZ*COS (BXAB*L) 

PM (19) =BXCD2*EIX*SINH(BXCD*L) -FZ*SINH (BXCD*L) 

PM (20) =BXCD2*EIX*COSH(BXCD*L) -FZ*COSH (BXCD*L) 

PM (9) =BYAB2*EIY*SIN (BYAB*L) +FZ*SIN (BYAB*L) 

PM (10) =BYAB2*EIY*COS (BYAB*L) +FZ*COS (BYAB*L) 

PM(ll) =-BYCD2 *EIY*SINH (BYCD*L) +FZ*SINH (BYCD*L) 

PM (12 ) =-BYCD2 *EIY*COSH (BYCD*L) +FZ*COSH (BYCD*L) 

PM(39)=-BP*EIP*COS (BP*L) 

PM (40) =BP*EIP*SIN (BP*L) 

7 CONTINUE 
RETURN 
END 

c 

SUBROUTINE QFORM ( W , L , Z , EIX , EIY , EAZ , EIP , MPL , IPL , FZ , AGKX , AGKY , QU , QS ) 
REAL L, MPL, IPL 
DIMENSION QU (40) , QS (40) 

W2=W*W 

IF (EIX-FZ*FZ* . 071) 1, 1, 2 

1 ARG=MPL/ (FZ-EIX*MPL*W2 /AGKX) 

BXAB=W*SQRT (ARG) 

BXCD=18 . 42 /L 
GO. TO 3 

2 CONTINUE 

AD= . 5* (MPL*W2 /AGKX-FZ/EIX) 

AE=MPL*W2/EIX 
ARG=AD+ SQRT ( AD *AD+AE ) 

BXAB=SQRT (ARG) 

ARG=-AD+SQRT (AD*AD+AE) 

BXCD= SQRT (ARG) 

3 CONTINUE 

IF (EIY-FZ*FZ* . 071)4,4,5 

4 ARG=MPL/ (FZ-EIY*MPL*W2 /AGKY) 



BYAB=W*SQRT (ARG ) 

BYCD=18 . 42 /L 
GO TO 6 

5 CONTINUE 

AD= . 5* (MPL*W2 / AGKY-FZ /EIY) 
AE=MPL*W2/EIY 
ARG=AD+ SQRT ( AD * AD+AE ) 
BYAB=SQRT (ARG) 

ARG=-AD+SQRT (AD*AD+AE) 
BYCD= SQRT (ARG) 

6 CONTINUE 
BZ=SQRT (MPL/EAZ ) *W 
BP=SQRT ( IPL/EIP) *W 
AD=1. 

C LINEAR DEFLECTION MATRIX 

CALL SET (QU, 3, 12) 

IF (Z) 8,9, 8 

9 QU ( 5 ) = 0 . 

QU (6) =1 . 

QU (7 ) =0 . 

QU (8) =1 . 

QU ( 2 1 ) = 0 . 

QU(22 ) =1 . 

QU (23 ) =0 . 

QU ( 2 4 ) = 1 . 

QU ( 3 7 ) = 0 . 

QU (38) =1 . 

GO TO 10 

8 CONTINUE 

QU ( 5 ) =SIN (BXAB*L) 

QU (6) =COS (BXAB*L) 

QU ( 7 ) =SINH (BXCD*L) 

QU ( 8 ) =COSH ( BXCD*L) 

QU (21) =SIN (BYAB*L) 

QU (22) =COS (BYAB*L) 

QU (23 ) =SINH (BYCD*L) 

QU (24) =COSH (BYCD*L) 

QU (37 ) =SIN (BZ*L) 

QU (38) =COS (BZ*L) 

10 CALL ADD (AD, QU, 0 . , QU, QU) 

C ANGULAR DEFLECTION MATRIX 

CALL SET(QS,3, 12) 

IF (Z) 11,12,11 
12 QS (17 ) =BXAB 

QS (18) =0 . 

QS (19 ) =BXCD 
QS (20 ) =0 . 

QS ( 9 ) =-BYAB 
QS (10 ) =0 . 

QS (11) =-BYCD 
QS (12) =0 . 

QS (39) =0 . 

QS (40) = -l . 

GO TO 13 

11 CONTINUE 
QS (17) =BXAB*COS (BXAB*L) 

QS (18) =-BXAB*SIN (BXAB*L) 

QS ( 19 ) =BXCD*COSH (BXCD*L) 

QS (20) =BXCD*SINH (BXCD*L) 

QS ( 9 ) = -BYAB * COS ( BYAB *L ) 

QS (10) =BYAB*SIN(BYAB*L) 

QS ( 11 ) =-BYCD*COSH (BYCD*L) 
QS ( 12 ) =-BYCD*SINH (BYCD*L) 
QS (39) =-SIN (BP*L) 

QS (40) =-COS (BP*L) 

CALL ADD(-1.,QS,0. ,QS,QS) 
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RETURN r 

END 

SUBROUTINE PQFORM (W, L, E, P, Q, DUM, DUN) 

REAL L, LI , MPL , IPL 

DIMENSION L(9),E(5),P(5),Q(5), DUM ( 5 ) , DUN ( 5 ) 

NBEAM=L ( 1 ) + . 01 

DO 1 IBEAM = 1 , NBEAM 

LI=L ( IBEAM+4 ) 

EIX=E( ( IBEAM- 1) *9+5) 

EIY=E( ( IBEAM- 1) *9+6) 

EAZ=E ( ( IBEAM- 1) *9+7) 

EIP=E ( (IBEAM-1) *9+8) 

MPL=E ( (IBEAM-1) *9+9) 

IPL=E ( (IBEAM-1) *9+10) 

FZ=E( (IBEAM-1) *9+11) 

AGKX=E ( (IBEAM-1) *9+12) 

AGKY=E ( (IBEAM-1) *9+13) 

CALL PFORM ( W , LI , 0 . , EIX , El Y , EAZ , El P , MPL , I PL , FZ , AGKX , AGKY , DUM 
1 , DUN) 

IF (IBEAM-1) 2, 3, 2 

3 CALL MAKE (P, DUM) 

GO TO 4 

2 CALL JUXTV ( P , DUM , P ) 

4 CALL JUXTV (P, DUN, P) 

CALL QFORM(W,LI, 0. , EIX, EIY, EAZ , EIP , MPL, IPL, FZ, AGKX, 

1 AGKY, DUM, DUN) 

IF (IBEAM-1) 5, 6,5 

6 CALL MAKE (Q, DUM) 

GO TO 7 

5 CALL JUXTV (Q, DUM, Q) 

7 CALL JUXTV ( Q , DUN , Q ) 

CALL PFORM (W, LI , LI , EIX, EIY, EAZ , EIP , MPL , IPL, FZ , AGKX, 

1 AGKY , DUM, DUN) 

CALL JUXTV (P, DUM, P) 

CALL JUXTV ( P , DUN , P ) 

CALL QFORM (W, LI, LI, EIX, EIY, EAZ , EIP , MPL, IPL, FZ , AGKX, 

1 AGKY, DUM, DUN) 

CALL JUXTV (Q, DUM, Q) 

1 CALL JUXTV (Q, DUN, Q) 

RETURN 

END 

c 

SUBROUTINE WSEARCH (W, DW, A, DETNEW, DETOLD, WI) 

DIMENSION A ( 5 ) , DETOLD ( 5 ) , DETNEW ( 5 ) , WI ( 5 ) 

707 FORMAT (7E12. 5) 

NA=A ( 1 ) 

A7 7 =A (77) 

CALL MAKE (DETOLD, DETNEW) 

DETO=l . 

DO 3 I =1 , NA 

3 DETO=DETO*DETOLD(I+4) 

CALL UPPER (A, DETNEW) 

DETN=1 . 

DO 1 1=1, NA 

1 DETN=DETN*DETNEW ( 1+ 4 ) 

FREQ=W/6 .2831853 
PRINT 707, FREQ , A7 7 , DETN 
IF (DETO*DETN) 4,2,2 

4 WROOT=W-DW-DETO*DW/ (DETN-DETO) 

NW=WI(1) +1.001 

WI (1)=NW 
WI (NW+4 ) =WROOT 
CALL SPIT (WI , 3H WI) 

2 RETURN 
END 



c 

SUBROUTINE UPPER (A, DETA) 

C GERNERATES THE DETERMINANT, DETA, OF MATRIX A 

DIMENSION A ( 5 ) , DETA ( 5 ) 

101 FORMAT (2E15. 7) 

N = A ( 1 ) + .01 
NM1 = N - 1 
NPl = N + 1 
DO 5 K = 1 , NMl 
KROW=K 

AKK=A(K*N+K-N+4) 

IF (AKK) 1,2,1 

2 KROW=KROW+ 1 

IF (KROW-N) 7, 7, 5 
7 AKK=A(KROW*N+K-N+4) 

IF (AKK) 3, 2, 3 

3 DO 4 J=K, N 
AD=A ( K*N+ J-N+4 ) 

A (K*N+J-N+4 ) =A (KROW*N+J-N+4 ) 

4 A (KROW*N+ J-N+4 ) =AD ' 

1 AKKl = 1 . /AKK 
Kl = K + 1 

DO 6 KI = K1,N 

AKIK = A (KI*N + K - N+4)*AKKl 
A (KI*N + K - N + 4) =0 . 

DO 6 L = Kl , N 

KIL = KI*N + L - N + 4 

A(KIL) = A (KIL) - AKIK*A (K*N + L - N + 4) 

6 CONTINUE 

5 CONTINUE 
AD=1 . 

CALL SET (DETA, N, 1) 

DO 16 I = 1 , N 

DETA ( 1+4 ) =A ( I*N + I - N + 4) 

EYE=I 

AD=AD*DETA ( 1+4 ) 

IF (AD) 26,27,26 
27 PRINT 101, AD, EYE 
AD=1 . 

2 6 CONTINUE 

16 CONTINUE 

RETURN 

END 

c 

SUBROUTINE DIAG ( A , SHAPE ) 

DIMENSION A (5 ), SHAPE (5) 

123 FORMAT ( 110, 5E12. 4) 

NA=A ( 1 ) + . 0 1 

NAM1=NA-1 

DO 1 I = 1 , NAMl 

IPl=NA-I+l 

AD=0 . 

DO 2 J = IP1 , NA 

2 AD=AD-SHAPE ( J+4 ) *A ( (NA-I-1 ) *NA+ J+4 ) 

SHAPE (NA-I+4 ) =AD/A ( (NA-I-1 ) *NA+NA-I+4 ) 

1 CONTINUE 

RETURN 
END 
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ABSTRACT 

During the past three decades the finite element method matured gradually and dominated 
almost exclusively all the engineering applications. But the practice of generating complex finite 
element dynamic models of aerospace structures has revealed a number of shortcomings. First, 
the high dimensionality of the models requires an order-reduction process before a control system 
can be designed. However, seemingly unimportant modes can be inadvertently eliminated which 
prove later to be significant to control system performance and stability. Second, use of finite 
element models for dynamic system identification is generally very difficult because of the 
extremely large amount of unknowns. Third, the structural damping is generally added ad hoc 
after generating an undamped model. Inaccurate damping results, especially for the modes which 
couple different types of motion. 

In contrast, distributed parameter models offer an alternative to finite element models for 
overall dynamic analysis and control synthesis for aerospace and aeronautical structures. First, 
the modal order does not have to be reduced prior to the inclusion of control system dynamics, 
which eliminates the risk involved with modal truncation. Second, distributed parameter models 
inherently involve fewer parameters, thereby enabling more accurate parameter estimation using 
experimental data. Third, it is possible to include the damping in the basic model, thereby 
increasing the accuracy of the structural damping. Recently, distributed parameter models have 
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been made for some large space structures. Distributed parameter models may provide an 
efficient complementary methodology to the finite element approach. 

The software package PDEMOD was initialized by the late Lawrence W. Taylor, Jr. at 
NASA Langley Research Center during the middle of the 1980’s. The initial interest in the 
package PDEMOD was to model the structural dynamics of general spacecraft configurations by 
using the distributed parameter approach, which consists of a three-dimensional network of 
flexible beams and rigid bodies. The building blocks from which three-dimensional configurations 
can be constructed consist of (1) beams, which have bending in two directions, torsion, and 
elongation degrees of freedom, and (2) rigid bodies, which are connected by any network of beam 
elements. The full six degrees of freedom are allowed at either end of the beam. Rigid bodies can 
be attached to the beam at any angle or body location. The modified Bemoulli-Euler beam 
equation is used to represent the bending, the wave equations for torsion and elongation. 

A system of partial differential equations (PDEs) is formulated and connected at the 
elements’ boundaries based on the compatibility conditions. The equations of motion for any 
number of rigid bodies are written in the frequency domain and in terms of the coefficients of the 
sinusoidal and hyperbolic functions which comprise the mode shapes. The force and moment 
vectors for both ends of a single beam element can be described in terms of the spatial derivatives 
of the solutions of the corresponding PDE’s. Distributed parameter models can therefore be 
generated for any three-dimensional configurations describable by PDEs joined at their 
boundaries. Because of Mr. Taylor’s sudden demise, it becomes an urgent task to summarize and 
sift his research achievements and make them available to the other researchers. We have 
continued his work and recovered the basic functions of this package which may provide an 
opportunity to more researchers to apply distributed parameter modeling techniques to a variety 
of aerospace structures. The verification of the code has been conducted by comparing the results 
with those examples for which the exact theoretical solutions can be obtained. 

By investigating the potential of the distributed parameter modeling technique, we expect 
that the PDEMOD may be further developed in the following aspects: (1) structural dynamics, 
modal frequencies and mode shapes; (2) parameter estimation of modal characteristics; (3) 
structural damping; (4) control system dynamics; and (5) design optimization. In its present 
stage, however, only the first of the functions has been completed and included in the package. 
Meanwhile, a massive effort is being conducted which is expected to modify the mathematical 
model and global system generating procedure to develop the methodology for control analysis 
purpose. Instead of using the coefficients of the solution functions, the transfer matrix may be 
used finally, which provides a much more convenient way to describe the state-vector transition 
from one point of the structure to the other. All these research outputs are planned to be contents 
of the modified version of the package. 


INTRODUCTION 

During the past three decades the finite element method had become extremely popular in 
a very wide field of engineering. As early as the 1960’s, the aircraft industry had developed in- 
house finite element programs. During the 1970’s, some general-purpose finite element programs 
such as NASTRAN were released for public use, bringing with them a significant technology base 
that led to development of numerous commercial finite element software systems. Later, the 
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various commercial packages were refined, and their technology base was expanded. NASA s 
research in computational structures technology (CST)[ij is helping to develop the finite element 
analysis to a new stage that is capable of using a general automated unstructured grid generation 
to discretize the aerodynamic field and the structure for thermal and structural analysis pj. 

But such a treatment having a very large number of elements is usually very expensive and 
time consuming. In static finite element analyses, models having over 10,000 degrees of freedom 
(DOF) are not uncommon. In fact, a finite element model of a helicopter fuselage currently under 
design contains about 27,000 grid points and 149,000 DOF. It is economically unfeasible and 
normally unnecessary to conduct dynamic analysis with this many unknowns. Most of the 
dynamic models based on the finite element approach must resort to modal truncation techniques 
to reduce the number of unknowns prior to dynamic analysis. However, the control spillover 
resulting from the modal truncation may lead to degradation of control system performance, or 
even create instability^]. 

The configuration of large space structures usually has the following common 
characteristics: extremely large dimension, light-weight design, high flexibility, rather uniform 
mass and stiffness distribution, low and closely spaced natural frequencies, and slight and/or 
improperly modeled damping. Structures with these characteristics are essentially distributed 
parameter systems by nature, and should be most accurately modeled as a continuously distributed 
mass and stiffness over the entire structural area rather than a sequence of finite mass elements 
coupled together as the finite element approach does. 

In contrast, distributed parameter modeling provides a very practical approach for overall 
dynamic analysis and control synthesis for aerospace and aeronautical structures. Recently, 
distributed parameter models have been made for some large space structures, such as, the 
Spacecraft Control Laboratory Experiment (SCOLE) [4] , Solar Array Flight Experiment (SAFE)pj, 
Space Station, Freedom (6] , Low-power Atmospheric Compensation Experiment (LACE) satellite 
model [ 7 ], and the Aerospace Large Flexible Manipulatory. The fundamental advantage of using 
the distributed parameter approach is to decrease the number of unknowns significantly. In the 
preliminary design stage, the detailed structural design is often neglected; therefore, a simple but 
efficient global structural model may be more beneficial to weigh the trade-off between the 
performance and cost, between the structural penalty and control system consummation, etc. 
Distributed parameter models may provide an efficient complementary methodology to the finite 
element approach. 

The software package PDEMOD was initialized by the late Lawrence W. Taylor, Jr. at 
NASA Langley Research Center during the middle of the 1980’s. The first release of his work on 
PDEMOD package was in 1987 (9 ,i 0 ,ii]. The initial interest in the package PDEMOD was to model 
the structural dynamics of general spacecraft configurations by using the distributed parameter 
approach, which consists of three-dimensional network of flexible beams and rigid bodies. The 
building blocks from which three-dimensional configurations can be constructed consist of (1) 
beams, which have bending in two directions, torsion, and elongation degrees of freedom, and (2) 
rigid bodies, which are connected by any network of beam elements. The full six degrees of 
freedom are allowed at either end of the beam. Rigid bodies can be attached to the beam at any 
angle or body location. The modified Bemoulli-Euler beam equation is used to represent the 
bending, the wave equations for torsion and elongation. 

A system of partial differential equations (PDEs) is formulated and connected at the 
elements’ boundaries based on the compatibility conditions. The equations of motion for any 
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number of rigid bodies are written in the frequency domain and in terms of the coefficients of the 
sinusoidal and hyperbolic functions which comprise the mode shapes. The force and moment 
vectors for both ends of a single beam element can be described in terms of the spatial derivatives 
of the solutions of the corresponding PDE’s. Distributed parameter models can therefore be 
generated for any three-dimensional configurations describable by PDEs joined at their 
boundaries. The manual labor of generating such models is therefore avoided. 

Because of Mr. Taylor’s sudden demise, it becomes an urgent task to summarize and sift 
his research achievements and make them available to the other researchers. We continued his 
work and recovered the basic functions of this package which may provide an opportunity to 
more researchers to apply distributed parameter modeling techniques to a variety of aerospace 
structures. The verification of the code has been conducted by comparing the results with those 
examples for which the exact theoretical solutions can be obtained, for instance, a simple beam 
with various boundary conditions, etc. 

By investigating the potential of the distributed parameter modeling technique, we expect 
that the PDEMOD may be further developed in the following aspects: (1) structural dynamics, 
modal frequencies and mode shapes; (2) parameter estimation of modal characteristics; (3) 
structural damping; (4) control system dynamics; and (5) design optimization. In its present 
stage, however, only the first of the functions has been completed and included in the package. 
Meanwhile, a massive effort is being conducted which is expected to modify the mathematical 
model and global system generating procedure[ 7 ,i 2 j to develop the methodology for control 
analysis purpose^,^. Instead of using the coefficients of the solution functions, the transfer 
matrix may be used finally, which provides a much more convenient way to describe the state- 
vector transition from one point of the structure to the other. All these research outputs are 
planned to be contents of the modified version of the package. 


PARTIAL DIFFERENTIAL EQUATIONS (PDEs) 

A complex large space structure can be decomposed into simple pieces. A network of 
distributed parameter elements and the attached rigid bodies are connected to represent the 
structural dynamics of the complex flexible spacecraft. Each flexible beam element exhibits lateral 
bending u and v in two axes, axial deformation w, and torsion \jr, as shown in Fig. 1, which can be 
independently described by a variety of PDEsfisj. 



Figure 1 A Beam Element 


The bending behavior of a beam element can be described by a modified Bemoulli-Euler 
beam equation which includes Euler bending stiffness, Timoshenko shear, and axial-force stiffness 
for lateral deflections in the x-z and y-z planes. The PDE for bending in the x-z plane is 
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( 1 ) 


mu + EI x u " + GAu" + F 0 u" = q x (z, t) 

The corresponding equation for bending in the y-z plane is 

mv + EI y v” + GAv' + F 0 v' = q y (z, t) (2) 

The axial and torsional dynamics are represented by wave equations. 


and 


mw - EAw” = F t (z, i) 


J Y y-GI y y' = M z (z,t) 


( 3 ) 

( 4 ) 


respectively. The PDEs provide the relationships between the modal frequency and the 
eigenvalues for the mode shape functions. The Euler and wave equations can be solved for the 
zero damping cases to produce the following relationships between the modal frequency to and 
the eigenvalue For bending in the x-z plane, 
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( 5 ) 


For bending in the y-z plane. 
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For elongation and torsion, we have 


and 


respectively. 




( 7 ) 

( 8 ) 


MODE SHAPE FUNCTIONS 

The solutions of these partial differential equations for zero damping produce the 
sinusoidal and hyperbolic spatial equations which comprise the mode shape functions. For the 
case that F 0 =0, the bending mode shape in the x-z plane is, 
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u = A x sin P x z + B x cos P x z + C x sinh P Z z + D x cosh p x z (9) 

Similarly, for the bending in the x-z plane the mode shape function is, 

v = A y sin p y z + B y cos P y z + C y sinh (5 y z + D y cosh (5 y z ( 1 0) 

The undamped mode shape function for elongation along the z-axis is 

w = A x sin P l z + B I cos P x z (11) 

The undamped mode shape function for torsion about the z-axis is 

y/= A y sin p y z +B y cos P Y z (12) 


These undamped mode shape functions are expected to be good approximations to the 
exact solutions for low level of damping. The mode shape of the entire configuration consists of 
all these functions, repeated for each beam element. Because of bending in two directions, 
elongation, and torsion, a total of 12 coefficients are needed for each beam element. A vector of 
the coefficients of these sinusoidal and hyperbolic functions is defined as the mode shape 
parameter vector, 

(0) = [A,B x C x D x A,B,C y D r A,B J A r B r f (1 3) 

The translational deflection vector is defined as, 


{(/} = M = [&(*)]{$} (14) 

w 

where, 

Qu Ql 2 Ql* Ql 4 o o o o o ooo 
[ 0 , 00 ] = 0 0 0 0 Ql 5 Ql 6 Ql 1 Ql* 0 0 0 0 

. 0 0 0 0 0 0 0 0 Ql 9 Ql’ ]0 0 0. 

the non-zero elements of the matrix [Q u (z)] are as follows, 

Ql 1 = sin fi x z Ql 2 = cos p x z Q™ = sinh P x z Q' u 4 = cosh p x z 

Ql 5 = sin P y z Ql 6 = cos P y z Ql 7 = sinhp y z Ql 8 = cosh j 3 y z 

Ql 9 = sin/3,z Ql 10 = cos p x z 

The angular deflection vector is defined as, 
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{£/'}= V ' = [&(z)]{0} (15) 

where, 

Ql 1 Ql 2 Ql 3 Ql 4 o o o oooo o 

[&(*)]= oooo a 25 o 26 a 27 a 28 ooo o 

.0 0 0 0 0 0 0 000 Q 3 /" Q 3 ' 3 _ 


the non-zero elements of the matrix [Q s (z)] are as follows, 


Ql' =(5 x cos p x z 
Qs 5 = Py COS (3 y Z 

Q]'" = sin/5^ 


Ql 2 =-(5 x sin P x z 
Qs & ~~Py SlTX Py Z 

Q)' n = COSPyZ 


Ql 3 = P x cosh p x z 
QT =P y cosh P y z 


Ql 4 = P x sinh p x z 
Q? = p y sinh p y z 


FORCES AND MOMENTS 

Next, it is necessary to express the forces and moments at either end of the beam element 
in terms of the mode shape parameter vector {0} . The force vector is as 

F x j \EI X W 

{F} = ’ F,' = * EI y A = [P f (z)]{0} (16) 

F t \ [EAw\ 

where, 

P" Pp 2 Pp Pp 0 0 0 0 0 0 0 0 

[P F (z)]= 0 0 0 0 Pf P 36 P 31 P 3 * 0 0 0 0 

.00000000 P 39 P 3 - ]0 0 0. 

the non-zero elements of the matrix [Pf(z)] are as follows, 

Pp = -PlEIx cos P x z P' 3 = P\EI X sin p x z 

Pp = P\EI X cosh p x z Pp 4 = P\EI X sinh p x z 

P 35 =~P\EI y cos P y z P 36 = P\EI y sin P y z 

Pp = P 3 y EI y cosh P y z P 3i = p 3 y E! y sinhP y z 

Pp 9 =P,EA cos P z z Pp'° = ~P Z EA sin p t z 

The moment vector can be expressed as 


7 



A/1 \EI x u' 

{M} = j M y ,► = EI y v” = [P M (z)]{e) ( 17 ) 

.A/J [g/,^ 

where, 

^ P^ 3 0 0 0 0 0 0 0 o' 

te/(*)] = 0 0 0 0 P m P m P m P m 0 0 0 0 

.0 0 0 0 0 0 0 000 Ph u ph u 

M M J 

the non-zero elements of the matrix [P M (z)] are as follows, 

P m = -PIET, sin fe P* = -fcEI x cos p x z 

P m = Pl EI * s inh p x z P" = p 2 x EI x cosh/3 x z 

P M = -P 2 yEIy sin p y z P* = -ft 2 EI y cospyZ 

P m = P] EI y sinh P y z P 2 J = p)EI y cosh p y z 

P M U = Py GI y cos PyZ P^ 2 = -p v GI v sin PyZ 

It is also necessary to account for changes in axes from each beam to the body to which it 
is attached, and for points of attachment at some distance from the center of gravity of the body. 
The forces and moments that the beam i applies to the bodyj expressed in the body j’s coordinate 
system are, 

iF},=m f {F\={T] t {P r \{e), (18) 

{M}„ =[r] ; ,({M),+M,{F),)=[r] J ,([^j (+ [/ ? ] j ,[/> F ]){«, (19) 

where, the coordinate-transformation matrix from the beam i’s coordinate system to the body j’s 
coordinate system is 

cos(Af,.,x,) coaiX^y,) cos 

[71,.,= cos (7,,x,) cos (Yj,y f ) cos (7,.,z,) (20) 

_cos(Z,,x,) cos(Z,,jy,) cos(Z,,r,)_ 

the eccentric matrix at the attachment point between the beam i and bodyj is 

0 -r z r r 

[P]„ = r z 0 -r x (21) 

~ r T r X 0 I 
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RIGID BODY MOTIONS 


A Newtonian or inertial frame of reference XoY 0 Zo is used for the motions of all beam 
elements and rigid bodies. 



Figure 2 Rigid Body Motion 


Consider body j connecting several beams, of which beam k is taken into account as a 
datum. From Fig. 2, we have 

at time t=t 0 : R a , = K, + h 

at time t=t: R a< = R a<j + u, = R Cq +r 0 + u t 

On the other hand, R = R +r t , so, 

R Cl = R„, ~r, =R Ct +r 0 +u t -r t = R c<> + u, +r 0 xu; (22) 

where, u t and u' are the translational and angular deflection vectors at the attachment point 
between the body j and the beam k, respectively, and the vector difference r Q - r, has been 
expressed as the vector cross product of r 0 and u', i.e. Ar =r 0 -r t =r 0 xu', which can also be 
written in matrix form as 

Ar x 0 —r z r r u' x 

At y * — /*2 0 -r x ' u' y " (23) 

.ArJ -r T r x 0 

Therefore, Eq.22 can be written in matrix form as 


9 



fc},, = ted,. 


(24) 


Differentiating Eq.24, we get the acceleration of the body j’s center of gravity (C.G.), 

{4J,=m*<ff>*+lKyri A {i<% (25) 

The angular acceleration of the body j is simply expressed as 

{£},= [7 (26) 


STRUCTURAL DYNAMIC EQUATIONS 

The equations of motion for the connected bodies and elements consist of blocks of terms, 
assembled in an order dictated by the body and beam indices. The mass times the acceleration of 
each body is related to the sum of forces caused by each beam element and each applied force. 

M y {i? co } = SlF},, +Yj{f} ]m +{#} (27) 

r m 

where, is the sum of the i-beam forces acting on the body j; X{/} ■„ * s the sum of the 

» m 

m-applied forces acting on the body j; {g} is the gravitational vector. It is similar for the moment 
equation, 

t n m 

where, is the sum of the i-beam moments acting on the body j; Xl 71 9 [7?] [ T] {F} ; , 

i i 

is the sum of the moments acting on the body j caused by the i-beam forces {F} ;I ; is 

n 

the sum of the n-applied moments acting on the body j; and ^,[T]^[R] jm [T) jm {f] m is the sum of 

m ** 

the moments caused by the m-applied forces acting on the body j. 

For the case of without applied forces and moments and neglecting gravity force, referring 
to Eqs.25 and 26, we derive the following two equations. From the force equation, 

N > ([7i A « A +i/ii A m,{sn,)=E^ (2») 

i 

From the moment equation, 

\J\M \{u% =S[{M} J , + [r],[2?]„(T],{F}„] (30) 

i 
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Let us express both sides of Eqs.29 and 30 in terms of the mode shape parameter vector, then, 
and 

-w I [y] / [7'], 1 [a] l {e} t =I((7']„[p„] i +[/)),[7'] J1 [i' F ]J(e}, ( 32 ) 

Eqs.31 and 32 are the structural dynamic equations for the most general configurations. To 
demonstrate the overall procedure more clearly, let’s assume that there is only one beam element 
attached to the body j, that is, i=k=l. In this case, Eqs.31 and 32 will be simplified as 

[4j{0} f -{O} (33-1) 

3x12 12x1 

[A u ] {0}, = {0} (33-1) 

3x12 12x1 

where, 

Two equations in Eq.33 may be combined as 

[A] (0), = (0) (34) 

6x12 12x1 

where the system matrix [A] consists of the two block matrices [Af] and [Am]- We can see that 
the number of equations is less than the number of unknown parameters. The difference is six. 
This is because there are six rigid-body degrees of freedom (d.o.f.s) for this particular one-body- 
one-beam system. We should have, therefore, six constraint equations to fix the six rigid-body 
d.o.fis. For most common cases, six boundary conditions provide six constraint equations, which 
can also be expressed in terms of the mode shape parameter vector {0}i, but the format of the 
equations depends on what the specific boundary conditions are. Superimposing the constraint 
equations into Eq.34, we obtain a new system dynamic equation 

l AIM, ={0} (35) 

12x12 12x1 

which has a full-rank system matrix [ A] for any values of circular frequencies except the natural 
frequencies. Based on the condition that Det\A\ = 0, the natural frequencies of the system can be 
determined. 

For general configurations, the structured dynamic equation may become very 
complicated, but the procedure to generate the equation is the same as stated. In general, for a 
structural system consisting of J bodies and I beams, the structural dynamic equation has the form 
of 
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( 36 ) 


U . I {«, = fo) 

(6*y)x(12*/)(i2./)xl 


The difference between the number of unknowns and the number of equations is exactly equal to 
the number of rigid-body d.o.f.s. To fix the rigid-body d.o.f.s, it is sometimes necessary to 
consider the compatibility conditions besides the boundary conditions. The following two special 
cases must be paid more attention. 

(1) For the case that a rigid body may have more than one beam element attached as shown in 
Fig.3 where two beams are attached to the body j, additional constraint equations must be added 
to the system equation, Eq.36, which appears as for this particular example in the figure, 

[A] {0} = {0} (37) 

12x24 24x1 

It is clear that 12 extra equations are needed to fit the two bodies’ rigid-body d.o.f.s. While the 
boundary conditions provides six equations, the other six equations must be found from the 
compatibility conditions. 



Figure 3 A Rigid Body Attached by Two Beam Elements 


To account for the continuity in translational deflection at the two attachment points, the 
constraint equation must be satisfied, 




(38) 


The constraint equation for ensuring continuity in the angular deflection at the two attachment 
points is 




(39) 


Expressing Eqs.38 and 39 in terms of the mode shape parameter vectors of the two beams, we 
have 


and 


(m„,[al„ -wwnJaijH, = (m,*[&] n -tfwrulaJjw),, 


(40) 

(41) 


Combining Eqs.40 and 41, we find the six additional equations. 
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(2) For the case that a beam connects two rigid bodies at its two ends as shown in Fig.4 
where beam i connects two bodies. More discussion must be addressed for this case. Apparently, 
it seems that the system equation, Eq.36, had a full-rank system matrix, since Eq.36 appears as for 
this particular example in the figure, 


U]{e}, = { 0} 

12X12 12x1 


body ji 


beam i 






1 body 


777777 


Figure 4 A Beam Element Connecting Two Rigid Bodies 


(42) 


But, it is wrong. The problem is that body jx and body j '2 connect to the beam i at different 
attachment points. Looking back to Eqs.14 to 17, and 31, we find that the components of the 
system matrix [A], such as [Q u (z)], [Q,(z)], [P F (z)], [P M (z)], are functions of the beam’s 
longitudinal coordinate z. Based upon the connection between body ji and beam i, a set of 
equations, which is the same as Eq.34, can be found, 

[A J} (z, = 0)] {0}, = {0} (43-1) 


Similar equation exists between body j 2 and beam i, 

[4 ; (r, = A-)]{9}, = {0} (43-2) 

6x12 12x1 

The two equations in Eq.43 are not, however, independent since the system matrices [Aji] and 
[Aj 2 ] are both related to the same beam element, beam i. It can be proven that these two matrices 
are related by a constant matrix [$] which can be derived from beam i’s PDEs, that is, 

[A(i,=I,)]=W[A(^0)] (44) 


We can only, therefore, choose one set of equations from Eq.43, and the other six equations must 
be fixed by the corresponding boundary conditions. 


VERIFICATION EXAMPLES 

1. EXAMPLE 1: Bending of a Cantilevered Beam. 
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Length of the Beam L=130.0in. 

Bending Stiffness EI=4xl0 7 Lb. in 2 . 

Mass per Length m=0. 09556 Lb.sec 2 /in. 


Theoretical formula for circular natural frequency: co n = a n^~jT = 1-210 6a n . The coefficients a„, 

and the theoretical natural frequencies and the corresponding results from PDEMOD are listed in 
Table 1. 


Table 1 Results of Example 1 


No. of Modes 

- an 

Natural Frequency 

Theoretical Value 

PDEMOD 

COn 

fn, HZ. 

fn/Hz. 

i 

3.52 

4.2613 

0.6782 

0.6775 

2 

22.0 

26.6332 

4.2388 

4.246 

3 

61.7 

74.6940 

11.8879 

11.890 

4 

121.0 

146.4826 

23.3134 

22.570 

5 

200.0 

242.1200 

38.5346 

38.490 


2. EXAMPLE 2: Bending of a Clamped-Clamped Beam. 


Length of the Beam 
Bending Stiffness 
Mass per Length 


L=130.0 in. 

EI=4xl0 7 Lb. in 2 . 
m=0. 09556 Lb.sec 2 /in. 


Theoretical formula for circular natural frequency: 


[eT 


1.2106a n . The coefficients a n , 


and the theoretical natural frequencies and the corresponding results from PDEMOD are listed in 
Table 2. 


Table 2 Results of Example 2 


No. of Modes 

a« 

Natural Frequency 

Theoretical Value 

PDEMOD 

COn 

f„, Hz. 

fn, HZ. 

i 

22.0 

26.6332 

4.2388 

4.31 

2 

61.7 

74.6940 

11.8879 

11.88 

3 

121.0 

146.4826 

23.3134 

23.29 

4 

200.0 

242.1200 

38.5346 

38.50 

5 

298.2 

361.0009 

57.4551 

57.50 


3. EXAMPLE 3: Bending of a Cantilevered Beam with a Tip-Mass M. 

Length of the Beam L=3.077 ft. 

Bending Stiffness EI=175.9644 Lb. ft 2 . 
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Mass per Length m— 0.037 Lb.sec 2 /ft (slug). 

Mass at the Tip M=0. 1 53 8 Lb.sec 2 /ft (slug). 

Equivalent Stiffness of the Beam k=3EI/L 3 =l 8. 1202 Lb/ft. 

Theoretical formula for the first circular natural frequency: 0), = 


i 


M+023m 


The theoretical 


natural frequency and the corresponding result from PDEMOD are listed in Table 3. 


Table 3 Results of Example 3 


No. of Modes 

Natural Frequency 

Theoretical Value 

PDEMOD 

(0„ 

fn, Hz. 

fn,Hz. 

i \ 

10.566 

1.682 

1.736 


4. EXAMPLE 4: Torsion of a Cantilevered Beam. 


Length of the Beam L=130.0 in. 

T orsional Stiffness GI y =4x 1 0 7 Lb . in 2 . 

Polar Moment of Inertia per Length jyL=0.2907 Lb. sec 2 . 

Theoretical formula for the circular natural frequency: co n 


GL 


-me. 


( JJL)l ? 


283.4745«. The 


theoretical natural frequency and the corresponding result from PDEMOD are listed in Table 4. 


Table 4 Results of Example 4 


No. of Modes 

Natural Frequency 

Theoretical Value 

PDEMOD 

(On 

f n , Hz. 

fn, HZ. 

i 

283.4745 

45.1164 

45.10 

2 

566.9490 

90.2328 

90.19 


CONCLUDING REMARKS 

The computer software package PDEMOD is being developed to model the structural 
dynamics of general spacecraft configurations by using the distributed parameter approach. This 
paper provides the detailed description of the theoretical background used in the PDEMOD 
formulation. A complex large space structure is considered as an assembly of flexible beam 
elements and rigid bodies. Each flexible beam element is represented by four independent partial 
differential equations which exhibit lateral bending in two axes, axial deformation, and torsion. A 
system of partial differential equations is then formulated and connected at the elements’ 
boundaries based on the compatibility conditions. The equations of motion for any number of 
rigid bodies are written in the frequency domain and in terms of the mode shape parameter 
coefficients. The deflections, forces and moments for both ends of a single beam element can be 
described in terms of the spatial derivatives of the solutions of the corresponding PDE’s, further 
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expressed in terms of the same set of mode shape parameter coefficients. Distributed parameter 
models can therefore be generated for any three-dimensional configurations describable by PDEs 
joined at their boundaries. The verification of the code has been conducted by comparing the 
results with those examples for which the exact theoretical solutions can be obtained. 

By investigating the potential of the distributed parameter modeling technique, we expect 
that the PDEMOD may be further developed in the following areas: (1) structural dynamics, 
modal frequencies and mode shapes; (2) parameter estimation of modal characteristics; (3) 
structural damping; (4) control system dynamics; and (5) design optimization. In present stage, 
however, only the first of these functions has been completed and included in the package. 
Meanwhile, a massive effort is being conducted which is expected to modify the mathematical 
model and global system generating procedure to develop methodology for the control analysis 
purpose. Instead of using the coefficients of the solution functions, the transfer matrix may finally 
be used, which provides a much more convenient way to describe the state-vector transition from 
one point of the structure to the other. All these research outputs are planned to be contents of 
the modified version of the package. It is also necessary to conduct additional testing to establish 
the accuracy of modeling different types of real space structures so as to move the approach from 
academic curiosity to a practical alternative for engineering design. 
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ABSTRACT 

The modeling, analysis and control of a small manipulator system are simplified by 
considering the links as rigid bodies. For large flexible manipulator systems, however, the 
flexibility of the links and the joint compliance must be considered. To describe the kinematics 
and the dynamic behavior of a flexible manipulator system, the common approach is to use 
Lagrange’s equations for both the rigid-body degrees of freedom (d.o.f.’s) and the dynamic 
deflection d.o.f.’s caused by the flexibility. The generalized coordinates are associated with the 
rigid-body d.o.f.’s of the links, and the modal coordinates associated with the flexibility d.o.f s. 
The consequence is that a set of highly-coupled and non-linear simultaneous partial differential 
equations is generated: These equations are so complex and lengthy that it is extremely difficult, 
if not impossible, to expand them manually even for a lower degree-of-freedom manipulator with 
a lower number of modes assumed. For the flexible manipulators with greater complexity, the 
dynamic analysis is literally forbidden by any practical manual symbolic derivations. The 
computer symbolic derivation of flexible manipulator dynamics was then suggested by several 
researchers. 

For simplifying the analytical process to a realistically acceptable extent, this paper 
conceives a new mathematical treatment for dynamic analysis of large flexible manipulator 
systems. The essence of the idea is to separate the kinematics and flexibility analyses as two 
independent but successive steps in a small time interval. Superposing the kinematic result and 
the flexibility effect, the summation is viewed as the initial conditions of the next instant motion, 
and the dynamic analysis succeeds to the next time interval. Repeating this process, the kinematic 
analysis accompanying flexibility effect is accomplished for a required time period. As usual, the 
kinematic analysis is based upon the rigid-body link assumption, and the Lagrange’s equations are 
set up for these less amount of rigid-body d.o.f.’s, which are manually manipulative. The 
flexibility analysis for certain configuration of the manipulator system is conducted by using the 
distributed parameter system approach, along with the application of the transfer matrix method. 
Since an extremely complex analytical chore is resolved into two relatively simpler problems, the 
complexity of the dynamic analysis of large flexible manipulator systems is mathematically 
simplified. To demonstrate the applicability of the proposed methodology, an end-effector 
vibration suppression problem for a large manipulator system has been investigated. The 
manipulator system studied in this paper is a similitude of a NASA manipulator testbed for the 
research of the berthing operation of the Space Shuttle to the Space Station. The computational 
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results show that the proposed method is very effective for end-effector vibration suppression of a 
large flexible manipulator system. 


1. INTRODUCTION 

The modeling, analysis and control of a small manipulator system are simplified by 
considering the links as rigid bodies^]. For large flexible manipulator systems, however, the 
flexibility of the links and the joint compliance must be considered. The limitations of the rigid 
link assumption in the formulation and analysis of large flexible manipulator dynamics were 
investigated extensively. Several formulations can be found in the robotics literature, such as, 
recursive or non-recursive Lagrangian assumed modep. 7 ], generalized Newton-Euler method^], 
and Lagrangian using Raleigh-Ritz method^], etc. To describe the kinematics and the dynamic 
behavior of a flexible manipulator system, the common approach is to use Lagrange’s equations 
for both the rigid-body degrees of freedom (d.o.f.’s) and the dynamic deflection d.o.f.’s caused by 
the flexibility. The generalized coordinates are associated with the rigid-body d.o.f.’s of the links, 
and the modal coordinates associated with the flexibility d.o.f s. The consequence is that a set of 
highly-coupled and non-linear simultaneous partial differential equations is generated. These 
equations are so complex and lengthy that it is extremely difficult, if not impossible, to expand 
them manually even for a lower degree-of-freedom manipulator with a lower number of modes 
assumed. For the flexible manipulators with greater complexity, the dynamic analysis is literally 
forbidden by any. practical manual symbolic, derivations. The computer symbolic derivation of 

• .flexible manipulator dynamics was later suggested by several researchers. Some of them wrote a 

• symbolic 1 manipulation program^], some of them suggested using the MATHEMATICA 
commercial software packagepoj. The basic functions of the symbolic manipulation program may 

.include symbolic simplification of. polynomials and rational expressions, linearization of 
trigonometric functions, automated evaluation of the relative significance of terms and neglecting 
the less significant terms, and even symbolic integration and differentiation. The application of 
computer symbolic derivation techniques alleviates the difficulty in the dynamic analysis of large 
flexible manipulator systems. 

For simplifying the analytical process to a realistically acceptable extent, this paper 
conceives a new mathematical treatment for dynamic analysis of large flexible manipulator 
systems. The essence of the idea is to separate the kinematics and flexibility analyses as two 
independent but successive steps in a small time interval. Superposing the kinematic result and 
the flexibility effect, the summation is viewed as the initial conditions of the next instant motion, 
and the dynamic analysis succeeds to the next time interval. Repeating this process, the kinematic 
analysis accompanying flexibility effect is accomplished for a required time period. As usual, the 
kinematic analysis is based upon the rigid-body link assumption, and the Lagrange’s equations are 
set up for these less amount of rigid-body d.o.f.’s, which are manually manipulative. The 
flexibility analysis for certain configuration of the manipulator system is conducted by using the 
distributed parameter system approach, along with the application of the transfer matrix 
methodpi]. Since an extremely complex analytical chore is resolved into two relatively simpler 
problems, the complexity of the dynamic analysis of large flexible manipulator systems is 
mathematically simplified. 

To demonstrate the applicability of the proposed methodology, an end-effector vibration 
suppression problem for a large manipulator system has been investigated. The manipulator 


2 



system studied in this paper is a similitude of a NASA manipulator testbed for the research of the 
berthing operation of the Space Shuttle to the Space Station^]. The system studied consists of 
two flexible links and three revolute joints which is assumed to be constrained in the vertical 
plane. There are only two rigid-body d.o.f.’s for this specific system. The two Lagrange’s 
equations are set up for the kinematic analysis, and the Runge-Kutta method is used to solve these 
non-linear partial differential equations numerically. The flexibility of the two links includes the 
lateral bending and axial elongation. The joint compliance is characterized by its torsional 
stiffness coefficient. The transfer matrices for the flexible arms and revolute joints have been 
constructed based on the partial differential equations. According to the compatibility conditions 
at the connecting points, the global system dynamic equation can be derived. From the 
corresponding boundary conditions, the characteristic equation for the global system is 
determined, from which the natural frequencies and mode shape functions can be found. 
Therefore, the transient response can be obtained. Joint moments are used as both displacement 
and control actuators. Control law computation proceeds in the frequency domain based on the 
pole-placement method^]. The computational results show that the proposed method is very 
effective for end-effector vibration suppression of a large flexible manipulator system. 

2. TWO-ARM FLEXIBLE MANIPULATOR SYSTEM 

The manipulator system (Figure 1) studied in this paper is a similitude of a NASA 
manipulator testbed for the research of the berthing operation of the space shuttle to the space 
station (Figure 2). This research testbed is planned to be the model of the berthing process 
constrained to move in the horizontal plane. Figure 2 illustrates the principal components of the 
facility. The Space Station Freedom (SSF) Mobility Base is an existing Marshall Space Flight 
Center (MSFC) Vehicle that has a mass of 2156.4kg. It represents a Space Station in the berthing 
operation. This vehicle is suspended on the MSFC flat floor facility using low flow-rate air 
bearings. The flexible appendage shown on the sketch will simulate solar panel disturbances. The 
vehicle has cold gas reaction jets to allow translational maneuvering. It also has a single gimbal 
for attitude control. The other vehicle, the Space Shuttle (SS) Mobility Base, is to be of similar 
construction and will be attached to the SSF Mobility Base with a flexible, two-link manipulator 
arm. The joints of the arms are driven by electric motors and are suspended by air bearings. 



Fig. 1 The Manipulator System Studied in the Paper 
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Fig.2 NASA MSFC Manipulator Testbed 


The system studied consists of two flexible links and three revolute joints: shoulder joint 
A, elbow joint B and wrist joint C. In the paper, it is assumed that the base frame XoY 0 Z 0 is fixed 
on the Shuttle assumed as a rigid body, with Xo-axis along the joint A axis. The orientation of the 
Y 0 and Z 0 axes about the joint axis Xo is chosen such that the resultant base frame forms a right- 
hand coordinate system. The frame XiyiZi for link 1 is defined as follows: the xi-axis coincides 
with the joint axis of the joint A, the Zi-axis is along the axis of link 1. The frame x 2 y 2 Z 2 is fixed at 
the joint B with x 2 -axis coinciding with the joint B axis, rotating with link 2, and z 2 -axis being 
along the axis of the link 2. Since the dimension of the end-effector can not be in general 
comparable with the dimensions of the two links, the end-effector will be abstracted as a rigid 
body represented by a mass point as a whole at the joint C. 

Each joint of the manipulator arms is driven by an individual actuator. The control 
moments x A , Tb, and Xc are also acting on the revolute joints A, B, and C, respectively. The joint 
compliance is characterized by its torsional stiffness coefficient. The corresponding input joint 
torques are transmitted through the arm linkage to the end-effector, where the resultant force and 
moment act upon the environment. The configuration of the corresponding rigid-body system of 
the one with flexibility can be specified by the two joint angles y/> and y/> as shown in Figure 1. 

3. KINEMATIC ANALYSIS UNDER RIGID-BODY LINK ASSUMPTION 


The common method for kinematic analysis under rigid-body link assumption in the 
robotics society is based on the solution of Lagrangian equation. 


dfdL y 
dt[dq t j 



(3.1) 


where, Lagrangian function L is defined as L-T-U, T and U are the kinetic and potential energies 
of the system, respectively. 0, is the generalized force corresponding to the generalized 
coordinate q t . For the two-arm system discussed in this paper, two rigid-body d.o.f.’s are chosen 
as the two angles % and y/ 2 as shown in Figure 1, where points D and E are the mass centers of 
the two arms. The kinetic and potential energies for both arms can be expressed as follows, for 
arm 1: 



/w i v o+^ / dV' i 


U x = ~™,gL x sin yr, 


(3.2) 
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For arm 2: 

T 2 =^m 2 vl+^I E yl = -m 2 [L]tf+L x L 2 y/,\j/ 2 cos(^, - y/ 2 ) + ~L\y/\] 

U 2= T W 2 g( 2L \ S ‘ n V\+ L 2 Sin Vi) 


(3-3) 


where, m„ L t are mass and length of the ith arm (i=l or 2). The Lagrangian function is then, 

L = T i +T 2 -U ] -U 2 

^ +\ m 2 m v\ + A A V\ v 2 cos (^i - ^2 ) + 3 A v\ ] 

+2w 2 )Z, sin^, +m 2 L 2 sin ^ 2 ] (3.4) 

Substituting Lagrangian function L into Eq.3.1 and taking corresponding derivatives, we derive 
the two Lagrangian equations: 


(^m, +m 2 )L\ij/ x +^m 2 L ] L 2 i// 2 cos(^, - y/ 2 ) + ^m 2 L } L 2 y/ 2 2 sin(^, - y/ 2 ) 
+K A y } -K B y/ 2 +|(m, +2 m 2 )gL, cos yr, = r, -r. 


(3.5-1) 


^-ot 2 AA<A cos(^, - yr 2 ) + -^w 2 AA^? sin (^i - ^ 2 ) 
+K B \y 2 +^m 2 gL 2 cosy/ 2 = t b -t c 


(3.5-2) 


Isolating ", and " 2 in each of the above two equations, Eq.3.5 can be reorganized as, 
[jot, +m 2 -\m 2 cos 2 ( j/, - r 2 )]A^l + T m 2 i i V'l sin2(^, - y 2 ) + \m 2 L x L 2 y\ sin(^, - y/ 2 ) 

31. cos(ii/ t - y/ 2 ) 1 3 

+ ^aV\ -[! + “ ~ 2 +( 2 «i +m 2)8 L i cos ^i -4 m 2^A cosj/ 2 cos(itq - Vi ) 


= (L< ~ x r)~ 


3L,(r 5 - r c )cos(^| - y/ 2 ) 


2L, 


(3.6-1) 


1 . % 2( 7 m,+/n 2 ) 

,^c«(r 1 -r2)-, 0M(ri . rt , 


+(ym ) + m 2 )Z^j ! tan(j/, + sin(yq -^ 2 ) 




1 + 


2 (^m ] +m 2 )L ] 

cos(^, -^) 


, (y/^+Wt^COS^ 

+ T( m I +2m 2 )gL, cosyq - 

r COS(^,-^ 2 ) 


= (*A ~ t b)~ 


2( T m, +m 2 )I | (T B -r c ) 
m 2 L 2 cos(|y, -i// 2 ) 


(3.6-2) 
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A set of two highly-coupled and non-linear simultaneous partial differential equations in Eq.3.6 
must be solved numerically. The well-known Runge-Kutta method is used. To do so, we define 


*1 = V\ (0. *2 = Vi (0, *3 = V 2 (0, *4 = V 2 (O' 

then, a set of four first-order simultaneous equations is generated, 

= /i (0 = *2 

x 2 = / 2 () = A x x\ sin2(x, -x 3 ) + 5,x 2 sin(x, -x^ + Qx, + Z),x 3 +E } cosx, 
+F X cosx 3 cos(x, - x 3 ) + G, cos(x, - x 3 ) + H x 

x 4 = / 4 (•) = A 2 x 2 tan(x, - x 3 ) + B 2 x] sin(x, - x 3 ) + C 2 x, + £> 2 x 3 + E 2 cosx, 


cosx, _ 1 

+E-, : r + G, — r + H-, 


where. 


__LnhA d _ T Shkh. 


A,=- 


2 cos(x,-x 3 ) 2 cos(x, -x 3 ) 2 

CQS(X, ~ X,) 


Kj 

A, 


1 7 | , C, =- , £,=-*-[ l+ '~ ‘ *' ] 


7 (m, +2m 2 )g£, 

£ 1 ” A ’ M “ 


*1 


2^2 


A, 


„ 3Z,,(r s r c ) ^ r, r B 

G i = “ •’ H i = ~ 




ILjA 


*i 


A, = I 3 /», +/n 2 _ T ot 2 cos (x, -x 3 )]i. 


and 


(3.7) 


(3.8) 


A=- 


(jw, +m 2 )Z, 2 


B 2 =- 


m 2 L x L 2 
2A-. 




G 2=-T^> ^2= , l 1 + 


K b 2( 3 m, +m 2 )I, 


"*2^2 


1 


e 2 =- 


(w, +2ffi 2 )gZ 1 ( 7 m, + m 2 )gL, 


2A, 




J 


1 

A 2 =[ 2 m 2 cos(x, -x 3 )-- -1I,A 2 

3cos(x, -x 3 ) 


^2 = ~ 


2( 3 w,+m 2 )Z 1 (r B -r c ) 


i/ 2 =- 


The iteration formula of the fourth-order Runge-Kutta method is as followS[i 4 ] : 


= *, j + 7 + 2* f 2 + 2k, 3 + 4 ) / = 1, 2, 3 and 4 

o 


where 


(3,9) 


^i. 1 “ ./i (*l.; > X 2.j > *3.; > *4,/ > 0 
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A, h h h h 

£(. 2 = fi( X \.J + ^k\ A ',X 2 J +~ki'), X 3 j +~^ 3 |,X 4 y + 2 ^ 4 . 1 >* + 2 ' 

h , /; /; h h 

kj,3 = + ^"^1.2> X 2 + ^k 2 2 , X 3j + 7^3,2 > X 4_; + 2 ^4.2 > * + 2^ 

^,4 = y"j (^1 j + ^\.3 > ^'2,/ ^^2,3 » ^3,/ ^^ 3,3 > *4,/ ^4,3 > ^ 

for the jth iteration, and h is the time interval. The accuracy of the method is in the order of A -5 . 
The motion of the rigid-body manipulator system can now be solved by using the iteration 
formula, Eq.3.9. 

4. FLEXIBILI TY ANALYSIS USING TRANSFER MATRIX METHOD 

In this paper, the inclusion of the flexibility of the manipulator arms is treated by the 
distributed parameter approach along with the transfer matrix method. The function of the 
transfer matrix is to relate the linear and angular deflections, forces and moments at one point in a 
structure to those at another point. The derivation of the transfer matrices for bending and 
elongation of a beam element is shown in Refs. 15 and 16, but is summarized here. The two 
flexible manipulator arms are represented by the Bemoulli-Euler equation and wave equation for 
their bending and elongation characteristics, respectively. The Bemoulli-Euler equation is in the 
form of 


&u v 1 d L u y 


dz 4 + a 2 y a 2 


= 0 


(4.1) 


where, a y 2 =k y /m, and k y =EI is the bending stiffness and m is the mass per length of the beam. The 
elongation vibration is described by the wave equation 


c fu . 1 cFu 


dz 1 


a 2 a 2 


- = 0 


(4.2) 


where, a z 2 =k z /m, and k z =EA is the axial stiffness. After separation of variables, Eqs.4. 1 and 4.2 
can be expressed in the spacial domain as 


ir r (z)-pp,(z)=o 

J/,(z) +#£/,« = 0 


(4.3) 

(4.4) 


where, P y 4 = o 2 /a y 2 , p z = co/a z , and co is the circular natural frequency. The solutions to Eqs.4. 3 and 
4.4 are 


and 


U y (z) = A sin P y z + B cos fi y z + C sinh P y z + D cosh P y z 
U z {z) = Ms\nP t z + NcosP t z 


(4.5) 

(4.6) 



where, A, B, C, D, M, N are the modal participation coefficients. If the state vector, (0), is 
defined, then 

{d) = [U y ,U„ ¥x ,F yy F t M x \ T (4.7) 

where, U y and U z are the displacements along y- and z-axes, y/ x is the rotary angle of the beam 
sections about x-axis, F y and F z are the shear and tensile forces respectively, M x is the bending 
moment about x-axis. The state vectors at the two ends of a beam element are related by a matrix 
[O], that is. 


{£(!)} = [O]{ 0 (O)} 


(4.8) 


where, the transfer matrix [<t>] involving bending and elongation of a beam element as described in 
Ref.16 is given in Eq.4.9, 


1 

— (COS PyL 
+cosh/?yZ.) 

0 


[4>] = 


1 


—Py(~ sin PyL 
+ sinh/7yL) 

k yP~y (SUl^?y L 


cos p z L 


2 % yyy\ 9Ul Hy* 

+ siiihfiyL) 

0 

1 2 

2 ^y^y ^ L 

+ cosh p y L) 


l 


Vy ^ p y L 

+ sinh^,I) 

0 

1 

— (COS PyL 
+ cosh PyL) 


1 


t(- smp v L 

2 k y py y 

+ sinh/?yZ,) 

0 

1 

i ( cosp v L 
2 kyfi* y 

+cosh PyL) 

1 


k yPy ( COS PyL (COS PyL 


-k 2 p 2 smp z L 


+co$hpyL) 

0 

1 

2 ^ yPy (” SATiPyL 
+ sinh/7yL) 


+ cosh/?^L) 
0 

+ sinh^yL) 


*zPz 


cos P Z L 


1 


si np z L 


y (-cos p v L 
VCyjTy 

+ cosh PyL) 

0 

+ sinh PyL) 

1 

~ Py (“ ^ 

+ sinhPyL) 

0 

1 

— (COS PyL 


+ cosh PyL) 

(4.9) 

in which the elements at the 2nd and 4th rows and columns reflect the elongation, the others for 
the bending. Because the two arms are not in the same orientation, it is necessary to account for 
the alignment of the two arms. The relationship of the two local coordinate systems can be 
described by a coordinate transformation matrix [T 2 i], that is. 


I y 2 


-fci- 


i* 


(4.10) 


where. 


[r„l = 


cosCv 2 ,^,) cos0v 2 ,z,) 
cos(r 2 ,y,) cos(z 2 ,Z|) 
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The transfer matrix of a revolute joint R has been derived in Ref. 11. A revolute joint R is 
abstracted as a massless torsional spring with spring constant kn> R through which two elements & 
and ‘b” are connected. The actuator fixed at the joint R will produce a control moment Tr. The 
state vector {0}, at the end of element Vis related to the state vector {9} b at the end of element 
“b” through the transfer matrix [<E>] of the joint R by 

< 41l > 

where, the transfer matrix [Or] of the joint R is given by 

kl 


1 

and the control-influence vector, 


(4.13) 




0 0 




0 0-1 



•Applying the general expression, Eq.4.1 1, to the joints A and B, we find out that, 


{0( Zl = 0)}„ = [O^HeUzr, = 0)}, 

(4.14) 

{«>, = £,)), = Mfe = 0)}, +M r . 

(4.15) 


where, the subscripts 0, 1 and 2 stand for the Shuttle Base, arm 1 and arm 2, respectively. 

5. SYSTEM DYNAMIC EQUATIONS FOR A SPECIFIC CONFIGURATION 

For a specific configuration at an arbitrary time instant, the system dynamic equation was 
derived in detail in Ref. 11, that is 

=£,)}, =U]{0f2, =0)l o +[B]{r! (5.1) 


where, the system matrix [A]=[ 02 ][ < I>b]" i [Oi][<I)a]’ 1 ; the control-influence matrix 

-{Be}]; the control vector {t}=[ t a , Tb, x c ] T 
The element transfer matrices [0*]’s and control influence vectors {B* } s are associated with the 
two links and the three joints, respectively, recognized by the corresponding subscripts. 

The dynamic properties at a certain configuration without control actions are the inherent 
properties of the system, which are varied with the change in configuration when the manipulator 
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system is in motion. The method to derive these inherent dynamic properties is straightforward. 
Consider the boundary conditions (B.C.’s) of the system in Figure 1, at the fixed end (attachment 
point to the Shuttle Base), 

u yt (r, = o) = U H (z, = 0) = ¥ Xo (z, = 0) = 0 (5.2) 

at the free end (end-effector), 

Fy, fe = A) = F h fe = L, ) = M XI (*, = L ,)-. = 0 (5.3) 

Applying the B.C.’s, Eqs.5.2 and 5.3, to Eq.5.1 without control action, two equations can be 
derived, 

; k(*,=°)| fcfe =■£,)' 

[4,1 F,(r,=o) =^,(^ = 4) (5.4) 

M,( z i=°)\ t kk=A)J, 

and 

>k= 0)' 

[■^22] ~ j = ® (5*5) 

A/,(z,=0)J o 

■ where, [A12] and [A22] are the block matrices of the matrix [A], The condition for Eq.5.5 having 
non-trivial solution is that 

DET[A 22 ]=0 (5.6) 

Eq.5.6 is the characteristic equation of the system from which the natural frequencies co's can be 
derived. After obtaining the natural frequencies from Eq.5.6, we can derive the mode shape 
functions from the equation below, 

[G]{Q=0 (5.7) 

The detailed derivation of Eq.5.7 can be found in Ref. 16. The vector {Q consists of the modal 
participation coefficients for the beams 1 and 2 appearing in Eqs.4.5 and 4.6, that is,{Q=[Ai, Bj, 
Ci, Di, Mi, Ni, A2, B 2 , C2, D 2 , M 2 , N2] t . Normalizing Eq.5.7 with N2=l, Eq.5.7 can be solved to 
obtain the modal participation coefficients for the two beams, thereby the mode shape functions 
can be obtained based on the solution equations, that is, Eqs.4.5 and 4.6, for the ith beam, 

U y , 0,) = A sin P y Zi +B, cos (3 y z t +C i sinh P y z t +£>, cosh^z,. (5.8) 

and U h (z,. ) = M t sin P 2> z, + N, cos& ( z, (5 . 9) 

It is assumed that the control actions are related to the feedbacks of the nodal displacements and 
velocities, that is. 


10 



r, = [*,]{«(*, =0)j, = [*,, + **S.** + V-** + VAO.0]{«U = 0)} ] (5.10) 

r, =[*,]{<(*, = 0)} 2 =[*„, +k,S,k, t +k B S,k, t +^J,0,0,oj{«(r 1 =0)} 2 (5.11) 

= [^ cl {^( Z 2 = ^2 )} 2 = [^ C , + *«£** + = A )}2 (^-^2) 

Inserting the control actions into Eq.5.1 and combining the similar terms, the closed-loop system 
dynamic equation follows, 

{^(z 2 =I 2 )} 2 =[I]{^z 1 = 0)} o (5.13) 


where, [^1 = [<!>,][<>,] '[*,][«>,] ‘.and 

[®„] =K]+M*„]. [®.l = K) + )[*.] . fa] = (w + fa] 

By applying the B.C.’s, Eqs.5.2 and 5.3, the closed-loop characteristic equation, DEl\ A 22 ] = 0, 
can be derived, where, [^ 2 ] is a block matrix of the matrix [A\, from which the closed-loop poles 
can be found. 

6. SUPERPOSING RIGID-BODY KINEMATICS AND FLEXIBILITY EFFECT 

This, paper conceives a new . mathematical treatment for dynamic analysis of large flexible 
manipulator systems. The essence of the idea is to separate the kinematics and flexibility analyses 
as two independent but successive steps in a small time interval. Section 3 of this paper gave the 
kinematic analysis assuming a rigid-body link based on the Lagrangian equation method, using the 
Runge-Kutta numerical approach. Sections 4 and 5 provided a method for system dynamic 
analysis due to flexibility at a specific configuration of the manipulator system. Compared with 
the macroscopic motion of the manipulator system, the motion resulted from flexibility is only a 
“microcosmic” motion. Only after a long-term effect is accumulated will the flexibility effect be 
significant. It allows us, therefore, to make the assumptions: at a certain configuration of the 
manipulator system, the deflections and the rates of deflection of the arms due to flexibility are 
small, and the elongation deformations are high-order infinitesimal so that they are neglected. 
Superposing the rigid-body motion and the motion due to flexibility, we have 

¥i( z i>0 = ¥i(0 + d ¥i( z u0 411(1 ¥2( z 2> t )='/'2( t ) +d V / i( z 2’ i ) ( 61 ) 

where, the small perturbation can be assumed as (cf. Fig.3), 

d V\( z \J) = —y\( z u t ) and = “T2( z 2>0 ( 6 - 2 ) 

z, z 2 

Therefore, the two generalized coordinates defining the instantaneous motions of the two arms in 
the Lagrangian equation, Eq.3.5, would be 
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A 


Fig. 3 Small Perturbation of the Angles 


y x (i) = \ { / x (z x ,t)-—y ] (z x ,t) and y/ 2 (t) = P 2 (* 2 .0 “ ^(* 2.0 ( 6 - 3 ) 

2, 2 2 

Substituting Eq.6.3 into Eq.3.6 and neglecting high-order infinitesimals, we can express the 
Lagrangian equations in terms of — , and ~ 2 as follows, 


= A \V\ sin2(^, -i// 2 ) + B x y/ 2 sin(^, -y/ 2 ) + C,^, + £, cos^, 

+F i cos \p 2 cos(^7, - ) + G, cos(^, - y7 2 ) + /f, + — .y, 


y/ 2 = A 2 y/ x tzn(y/ x -y/ 2 ) + B 2 y/ 2 2 sm(y/ x -y/ 2 ) + C 2 y/ x +D 2 y/ 2 +E 2 cos 

cosu7, 1 1 

+-^2 7= =“7 + G 2 tzt — N + y 2 

cos(((/| - V/,) COS( (Z, - (/,) Zj 


(6.4-1) 


(6.4-2) 


Defining w x = y/ x , w 2 = y/ x , w 3 = y/ 2 , and w 4 = y/ 2 , a set of four first-order simultaneous 
equations is obtained which is similar to Eq.3.8, 


w, = w 2 

w 2 = sin2(w, -w 3 ) + 5,w 4 sin(w, - w 3 ) + C,m', + D x w 3 + E x cosw, 

+F X cos>v 3 cos(m' ) - w> 3 ) + G, cos^, - w 3 ) + H x + —y x 

z i 


h< 4 = y4 2 w 2 tan(w, -■w' 3 ) + .6 2 m' 4 sin(w, -w 3 ) + C 2 w, + D 2 w 3 + ^2 cosw, 


COSM', _ 

+F 2 ZFJF. — + G 2 


cos^, - w 3 ) cos(w, - w 3 ) 


+H 2 +—y 2 


(6.5) 


The only difference between Eq.6.5 and Eq.3.8 is the inclusion of y x and y 2 , which represents the 
effect of flexibility, and can be solved based on the formulation described in Sections 4 and 5 as 
long as the instantaneous configuration is specified and the motion at the end of previous time 
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interval is known. By using the Runge-Kutta procedure in Section 3, Eq.6.5 can be solved 
numerically. A solution is thus obtained which represents the superposition of rigid-body 
kinematics and flexibility effect. 


7. SIMULATION EXAMPLES 

As mentioned earlier, the manipulator system studied here is a similitude of a NASA 
MSFC manipulator testbed for the research of the berthing operation of the Space Shuttle to the 
Space Station. The mechanical properties of the system are listed in Table 1. 

Table 1 Mechanical Properties of the Two-arm Manipulator System 


Beam Length, L, (in. 


Sectional Area, A, (in. 


Second Moment of the Cross Section, I, (in. 4 ) 


Modulus of Elasticity, E, 


Mass per length, (slug/in. 


Beam 1 

Beam 2 

60.0 

120.0 

50.27 

12.57 

107.0 

107.0 

201.06 

12.57 

0.1526 

0.0382 



The initial conditions of the system are assumedat y/ x = 90°, y/ 2 = 0°, y/ x =0, and y/ 2 = 0 . 
First, the difference between the results of under rigid-body link assumption and with flexibility 
effects are noted by three examples as shown in Fig.4. These examples exhibit the motions of the 
manipulator system from the specified initial conditions under the actions of gravity force, inertia 
force, and the constant joint control moments t a , t b , and x c as shown in each figure. The solid 
lines represent the motions under rigid-body link assumption, while the dashed lines represent the 
motions with flexibility effects. For clarity, the vibratory wave shapes for flexible beams were 
neglected in the figures. 

Next, several examples demonstrate the effectiveness for end-effector vibration 
suppression. The motion of the manipulator system is a continuous process, which stimulates 
vibration of the system at every time instant. The vibration suppression action continues 
throughout the whole process without interruption. The joint control moments will always 
change themselves based upon the control law given in Eqs.5.10 to 5.12. They play the roles of 
both actuating the manipulator system to fulfill a certain task and alleviating vibratory fluctuation 



13 
















Fig.4 The Difference Between the Results Of Under Rigid-Body Link Assumption 

And With Flexibility Effects 

Figs. 5 to 7 demonstrate the effectiveness for the end-effector vibration suppression at 
instantaneous positions 1, 2, 3 as shown in Fig.4, with the initial joint control moments Ta = 1 ft- 
klb, Xb= 2 ft-klb, and Xc=4 ft-klb. Both time histories without control (left) and with control 
(right) are shown in the figures for comparison. The vibratory motion of the end-effector is 
described in the second link’s coordinate system, x 2 y 2 Z 2 , as defined in Section 2. The upper two 
figures in Figs. 5 to 7 represent the results in y 2 -direction, the lower two in z 2 -direction. The 
instantaneous vibration can be suppressed in about 0.3 second for all positions studied. 

8. CONCLUDING REMARKS 

A new mathematical treatment for dynamic analysis of large flexible manipulator systems 
is derived. An extremely complex analytical chore is resolved into two relatively simpler 
problems, the complexity of the dynamic analysis of large flexible manipulator systems is, 
therefore, mathematically simplified to a realistically acceptable extent for practical manual 
symbolic derivation of the equations of motion. Since the equation of motion is a set of highly- 
coupled and non-linear simultaneous partial differential equations, the Runge-Kutta numerical 
procedure has been used to solve the equations. As an example, the vibration suppression 
problem of a similitude of a NASA MSFC manipulator testbed has been investigated in the paper. 
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The computational results show that the proposed method is very effective for end-effector 
vibration suppression for a large flexible manipulator system. 
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Fig.5 End-Effector Vibration Suppression at Instantaneous Position: \j/\ = 103.2°, y/ 2 = 3.6° 
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Fig.6 End-Effector Vibration Suppression at Instantaneous Position: y/\ = 121.8° , y/ 2 = 10.2° 
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